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Long-term  objective  of  our  research  effort:  A  systematic  hierarchical  design  approach 
for  advanced  rocket  propellants  that  exploits  intrinsic  and  engineered  anisotropy 

The  efficacy  of  new  propellants  based  on  hierarchically  anisotropic  materials  depends 
on  structures  in  the  unreacted  material,  but  also  on  the  influence  of  these  structures  in  the 
decomposing  and  burning  propellant.  Therefore,  we  think  it  is  useful  to  think  of 
anisotropy  at  three  levels: 

1)  Intrinsic  anisotropy  at  the  molecular  up  to  the  continuum  microscale  for  pure 
constituents 

2)  Manufactured  nano-  and  microscale  anisotropy  that  is  induced  by  manufacture 
specifications  of  the  composition  that  identifies  the  material  selection,  gross 
stoichiometry,  morphology,  and  constituent  placement  and  orientation 

3)  Mesoscale  anisotropy  persistence  during  the  physico-chemical  structural 
decomposition,  mixing,  and  reactive  processes,  templated  by  item  2),  that  results 
in  locally  anisotropic  effects  that  appear  homogenized  at  the  engineering  device¬ 
level  scale  that  ultimately  defines  the  performance  characteristics  of  the 
propellant. 

We  have  longstanding  interests  in  the  fundamental  physics,  chemistry,  and 
engineering  required  to  realize  the  anisotropy-  and  hierarchical-design-based  approach 
outlined  just  above.  The  goal  of  our  tri- institutional  project  is  a  focused,  integrated 
program  of  fundamental  theoretical  and  computational  research  designed  to  yield 
understanding  and  predictive  capability  regarding  anisotropic  thermal  transport  and 
energy  release  in  advanced  rocket  propellants.  While  fully  realizing  this  goal  is  not 
realistically  achievable  in  a  three-year  effort,  the  research  we  are  performing  will 
contribute  to  a  longer-term  practical  capability  for  the  a  priori  design  of  smart,  functional 
propellant  materials  based  on  intelligently  designed  polymer  nanocomposite  formulations 
augmented  by  non-traditional  additives  or  passivation  agents,  possibly  coupled  to 
external  fields. 

We  believe  that  the  optimal  approach  is  a  hierarchical  theoretical  framework  that 
combines,  with  as  much  rigor  as  possible: 

•  Fundamental  information  from  atomic-scale  simulations 

o  Nanoparticles  and  bulk 
o  Interfaces 

•  Detailed  continuum-based  mesoscopic 
models  of  interfaces  and  interphases, 
phase  transformations,  multi-phase 
transport,  and  chemistry 

•  Explicit  microstructure-resolved  RVE 
scale  and  larger  simulations  that 
generate  homogenized  results  induced 
by  mesoscale  constituents 


'lAto~250nm  ~100  nm  to  ~100  n  ~100nto~cm 


Schematic  depiction  of  nominal  spatial  scales, 
physical  situations  and  phenomena,  and  eritieal 
intraseale  and  seale-bridging  eonneetions  to  be 
studied. 
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•  Development  of  approaehes  that  define  engineering-scale  models  for  propellant 
combustion  uniquely  defined  by  the  hierarchical  structures  from  the  molecular, 
micro-,  and  mesoscopic  models.  These  engineering- scale  models  are  the  ones  that 
can  be  used  to  make  system  performance  estimates. 

Atomic-scale  studies  -  molecular  dynamics  (MD)  and  electronic  structure  theory  - 
are  providing  critical  knowledge  of  the  details  of  the  physics  and  chemistry  of  materials 
that  occur  at  and  in  the  vicinity  of  interfaces  in  advanced  propellant  formulations.  Our 
studies  are  resulting  in  consistent  identification,  characterization,  and  quantification  of 
the  unit  processes  of  anisotropic  mechanical  response,  transport  phenomena,  phase 
stability,  and  chemical  transformation  within  well-defined  (and  computationally 
tractable)  physical  situations. 

We  are  developing  mesoscopic  continuum-based  models  for  chemically  reactive 
materials.  The  general  approach  taken  follows  standard  methodologies  used  in  the 
combustion  community,  but  include  the  complexities  associated  with  anisotropy  and 
heterogeneity  that  must  be  considered  to  describe  nano-  and  microscale  physical 
processes  and  their  effects  on  chemistry.  Those  include  consideration  of  phase  changes 
including  solid-state  polymorphism  and  transitions  among  solid,  liquid,  and  gas;  and 
tracking  of  key  species  diagnostic  of  or  critical  to  the  overall  reaction  kinetics. 
Multiphase  diffusion  of  species  and  diffusion  of  oxidizing  and  metal  ions  in  oxide  melt 
solutions  (often  thought  to  be  involved  in  metal,  intermetallic,  and  thermitic  reactions)  is 
being  included.  The  mesoscopic  models  predict: 

•  Multiphase  decomposition  and  reaction  during  chemical  energy  release  at  initially 
separate  diffusive/reactive  interfaces 

•  Volumetric  reaction  and  gaseous  combustion 

•  Formalism  and  constitutive  theory  used  to  describe  reactants  and  their  products 
that  reflect  the  underlying  material  mechanics  and  physical  chemistry  of  the 
constituents 

We  are  developing  RVE-scale  methods  for  simulating  experimentally  relevant 
configurations  such  as: 

•  1-,  2- and  3-D  laminates  and  plies 

•  Both  periodic  and  stochastic  microstructure  with  designated  distributions  of 
materials. 

The  goal  is  to  simulate  a  piece  of  material  sufficiently  large  to  yield,  in  a  proper 
homogenization  over  the  "microscopic"  ensemble,  statistical  distributions  of  properties 
that  can  be  used  in  macro-scale  models. 

Our  research  is  primarily  fundamental  with  the  aim  to  provide  guidance  for  practical 
applications.  We  are  focused  on  establishing  tight,  direct  coupling  among  the  five  main 
elements  of  the  project  (see  the  figure  on  the  preceding  page).  The  fundamental  methods 
and  analysis  tools  we  are  developing  should  be  general,  so  we  do  not  anticipate  that  they 
will  be  successful  for  one  category  of  material  but  fail  for  another. 
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University  of  Missouri-Columbia  (MU) 

The  goals  of  the  MU  sub-project  are  to  (1)  develop  methods  and  models  for  atomic- 
scale  simulations  of  highly  anisotropic  propellant  constituent  materials  and  (2)  apply 
those  methods  to  the  careful  prediction  of  the  pressure-volume-temperature  equation  of 
state,  pressure-  and  temperature-dependent  crystal  and  liquid  thermal  and  transport 
properties,  and  fundamental  crystal  inelastic  deformation  behaviors  (plasticity  in  the 
basal  plane  and  the  nanoindentation  response  for  oriented  crystals).  We  also  studied  the 
melting  process  and  calculated  the  melting  curve  Tme\t=Tme\t{P)  for  a  range  of  propellant¬ 
relevant  pressures.  To  accomplish  these  goals  we  used  a  combination  of  molecular 
dynamics  (MD)  and  molecular  mechanics.  The  focus  of  our  effort  is  prediction  of 
properties  that  can  be  used  directly  in  the  larger  scale  theoretical  models  under 
development  at  the  University  of  Illinois  at  Urbana-Champaign  and  the  California 
Institute  of  Technology,  and  indeed  practically  every  property  we  discuss  below  was 
chosen  for  study  based  on  discussions  with  those  collaborators  (and,  in  the  second  year  of 
this  project,  has  been  used  by  them  -  Prof  Stewart  at  UIUC  in  particular). 

MU  First  Year  Summary  Comments:  Although  the  energetic  crystal  1,3,5-triamino- 
2,4,6-trinitrobenzene  (TATB)  is  not  of  practical  use  as  a  propellant  ingredient,  we  began 
our  work  with  this  compound  as  it  crystallizes  into  the  lowest  symmetry  crystal  class 
(triclinic)  and  exhibits  the  greatest  structural,  mechanical,  and  thermal  anisotropy  of  any 
energetic  compound  of  which  we  are  aware.  Thus,  assuming  access  to  a  suitable  potential 
energy  function  with  which  to  perform  the  simulations,  the  methods  we  are  developing 
should  be  applicable  to  practically  any  other  energetic  compound  of  interest.  We 
emphasize  that  our  studies  under  this  project  are  focused  on  the  development  of  an 
internally  consistent  thermo-mechanical  description  of  TATB  under  “equilibrium” 
conditions  rather  than  shock  loading. 

We  also  performed  high-level  (essentially  benchmark-quality)  electronic  structure 
(a.k.a.  “quantum  chemistry”)  calculations  of  proton  exchange  and  transfer  in  ammonium 
nitrate  (AN)  and  ammonium  perchlorate  (AP),  and  calculations  of  pathways  and 
energetics  for  the  overall  reaction  A1  -i-  CO2  ^  AlO  +  CO.  Proton  transfer  is  often 
thought  to  be  an  early  reaction  in  the  decomposition  of  inorganic  salts  such  as  AN  and 
AP,  which  are  commonly  used  oxidizers.  As  an  initial  stage  of  developing  chemical 
decomposition  schemes  for  AN  and  AP,  calculations  were  performed  to  study  the 
transformation  from  the  stable  acid-base  pair  for  isolated  formula  units  to  stable  ion  pairs. 
In  addition  to  calculations  for  pure  (AN)n  and  (AP)n,  the  effects  of  NH3,  AIH3,  and  BH3 
molecules  on  cluster  acid-base/ion-pair  stability  were  studied.  We  became  interested  in 
the  AI-1-CO2  reaction  chemistry  both  due  to  its  relevance  in  combustion  and  because  the 
combined  experimental  and  theoretical  descriptions  of  low-energy  AICO2  complexes 
appeared  to  be  inconsistent  and  incomplete.  Our  work  on  AI-1-CO2  should  provide  a 
complete,  and  consistent,  description  of  the  various  reaction  pathways  and  energetics  at 
the  fully  optimized  CCSD(T)  level  of  theory. 

More  explicitly,  the  properties  that  have  been  studied,  this  year  and  in  the  first  year  of 
the  project  (which  was  separated  from  the  current  funding  increment  by  a  gap  of  several 
months),  include: 

•  Force  field  development  for  TATB  (necessary  for  MD  simulations) 
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•  Anisotropic  thermal  conductivity  in  triclinic  crystalline  TATB  as  functions  of 
temperature  and  pressure 

•  TATB  liquid  state  transport  coefficients — thermal  conductivity,  shear  viscosity, 
and  self-diffusion — as  functions  of  temperature  and  density 

•  Second-order  elastic  tensor  for  TATB  for  finite  temperatures  and  pressures 

•  Calculations  of  plasticity  in  the  TATB  basal  plane  {i.e.,  in  the  (001)  plane) 

•  Calculations  of  nanoindentation  of  TATB  normal  to  the  (100),  (010),  and  (001) 
crystal  surfaces 

•  Melting  point  of  TATB,  Tmeit,  both  at  normal  and  elevated  pressure  {i.e.,  the  melt 
curve  rmelt=Tmelt(T’))  to  2.0  GPa 

•  Anisotropic  energy  transfer  from  one-dimensional  “hot  spots”  in  oriented  TATB 
crystals,  including  a  study  of  the  time  and  space  scales  on  which  energy  transport 
at  the  nanoscale  can  be  accurately  described  using  the  diffusive  heat  equation 
{i.e.,  the  continuum  heat  transfer  equation) 

•  Electronic  structure  predictions  of  the  transformation  from  acid-base  pairs  {e.g., 
nitric  acid  and  ammonia)  to  ion  pairs  {e.g.,  NH4^  and  NO3'),  that  is,  proton 
transfer,  in  small  clusters  of  ammonium  nitrate  (AN)n  and  ammonium  perchlorate 
(AP)n  containing  n  >  1  formula  units  of  AN  or  AP 

•  Electronic  structure  calculations  of  the  oxidation  of  A1  by  CO2. 

•  We  also  developed  and  published  a  method  by  which  the  data  storage 
requirements  for  saving  MD  simulation  trajectories  can  be  drastically  reduced,  by 
storing  data  in  a  stepwise-uniform  fashion  in  the  log(t)  domain,  in  a  way  that 
remains  consistent  with  the  Sharmon-Nyquist  sampling  theorem. 

MU  Second  Year  Summary  Comments:  In  the  second  year  of  the  project  we  completed 
research  that  carried  over  from  the  first  year  (and  from  the  “original”  first  year  that  was 
funded  under  a  different  but  intellectually  related  grant).  We  also  undertook  simulations 
of  the  mechanical  and  transport  properties  for  HMX  and  RDX.  These  are  both  propellant¬ 
relevant  materials,  and  indeed  HMX  is  a  monopropellant.  HMX  was  the  focus  of  the 
UlUC  effort  in  the  second  year  of  the  project,  and  many  of  the  properties  that  we 
obtained  for  HMX  are  presently  in  use  in  the  UlUC  simulations.  More  explicitly,  the 
tasks  completed  and  properties  that  have  been  studied  at  MU  in  the  second  year  include: 

•  Completion  of  manuscript  preparation  or  publication  for  several  of  the  items 
listed  above.  The  MU  reference  list  (below)  has  been  modified  to  provide  final 
citations  where  applicable.  Papers  that  were  submitted  in  the  first  year  but 
appeared  in  print  during  the  second  year  are  denoted  by  an  asterisk.  Papers  that 
were  listed  as  in  preparation  at  the  time  of  the  first-year  report  but  were  submitted 
or  published  in  the  second  year  are  treated  as  second-year  publications  and  thus 
not  denoted  by  an  asterisk. 

•  Calculation  of  the  complete  second-order  elastic  coefficients  for  P-HMX  on  the 
cold  curve  for  pressures  less  than  30  GPa.  (Manuscript  in  preparation.) 

•  Nanoindentation  experimentally  important  planes  in  P-HMX,  (110)  and  (001). 
(Work  in  progress  as  of  this  writing.) 
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•  Thermal  conduction  along  the  experimentally  relevant  p-HMX  crystal  directions 
(specifically,  directions  normal  to  the  experimentally  important  crystal  surfaces 
(110),  (Oil),  and  (010)). 

•  Simulations  of  Kaptiza  resistance  (i.e.,  thermal  interface  resistance)  of  various 
energetic  crystal  binary  interfaces  including: 

o  (001)1(100)  in  TATB  (of  interest  due  to  the  very  high  structural 
anisotropy) 

o  (01 1)1(1 10)  in  p-HMX  (a  grain  boundary  between  the  two  most  prominent 
interfaces  in  P-HMX) 

o  An  a-RDX/p-HMX  heterophase  interface 

o  A  nitromethane  (100)/(001)  grain  boundary 

o  All  of  these  were  in  progress  at  the  conclusion  of  the  project. 

•  Generalized  definitions  for  specifying  geometries  and  calculating  grain  boundary 
energetics  for  simple  (i.e.,  pure  twist)  and  more  complicated  grain  boundaries  and 
heterophase  interfaces.  (In  progress  at  the  end  of  the  project.) 

Note:  The  nanoindentation,  thermal  conduction,  and  Kapitza  resistance  studies  listed 
just  above  were  only  possible  the  benefit  of  a  method  published  by  us  very  recently 
{The  Generalized  Crystal  Cutting  Method  (GCCM),  Kroonblawd  et  ah,  Computer 
Physics  Communications  207,  232  (2016),  supported  by  a  different  AFOSR  grant] 
that  enables  the  construction  of  properly  periodic  simulation  cells  for  arbitrarily 
oriented  crystal  surfaces  or  oriented  homophase  or  heterophase  interfaces. 

Most  of  the  efforts  mentioned  above  are  described  in  somewhat  greater  detail  in  the 
following  paragraphs. 

I.  FIRST  YEAR  WORK,  SOME  OF  WHICH  WAS  COMPLETED  IN  YEAR  2 

■  Development  of  a  flexible  molecule  force  field  for  TATB:  We  developed' 
potential  terms  to  extend  an  existing^  all-atom  l,3,5-triamino-2,4,6- 
trinitrobenzene  (TATB)  force  field  (hereafter,  TATB  FF)  to  approximately 
reproduce  experimental  vibrational  spectra  and  normal  modes  determined  from  a 
density  functional  theory  calculation.  The  triclinic,  layered  crystal  packing 
structure  for  TATB  is  shown  in  Fig.  1.^  The  TATB  FF  has  been  validated  against 
experimental  determinations  of  the  temperature  and  pressure  dependent  lattice 
parameters,  heat  of  sublimation,  and  bulk  modulus.  ’  ’  Many  other  properties  of 
TATB  crystal  have  never  been  determined  experimentally,  due  in  large  part  to  the 
difficulties  in  growing  large,  high-quality  single  crystals.  The  extended  force  field 
is  currently  the  best  available  and  most-studied  flexible  molecule  model  for 
TATB  and  has  been  used  to  make  numerous  predictions  for  anisotropic 
thermomechanical  properties,  ’  ’  ‘  the  phase  diagram,  ’  transport 
coefficients, and  for  up-scaling  to  parameterize  engineering-scale  TATB 
models.^ 

■  An  improved  force  field  for  high-temperature  simulations:  A  repulsive 
intramolecular  potential  was  developed"'  to  remedy  inaccuracies  in  the  radial 
distribution  function  for  TATB  at  temperatures  (J)  above  750  K  predicted  by  the 
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TATB  FF.  This  potential  was  fit  to  data  obtained  from  eonstrained  geometry 
optimizations  in  GAUSSIAN  09  at  the  M06-2X/6-311G**  level.  The  improved 
force  field  predicts  realistic  intermolecular  0-H  distances  at  r>  750  K. 

Reducing  data  storage  costs:  Strategies  for  non-uniform  sampling  of  simulated 
relaxation  phenomena  were  developed'*^  that  can  reduce  data  storage  costs  by  over 
90%. 


Elastic  anisotropy  and  plastic  deformation  mechanisms:  Complete  elastic 
stiffness  tensor  of  TATB  was  determined^  at  pressures  (P)  of  1  atm  and  5  GPa  by 
fitting  second-degree  polynomials  to  energy-density-strain  curves  obtained  from 
molecular  statics  (MS)  simulations.  Pronounced  anisotropy  in  the  elastic  response 
is  predicted  with  significant  stiffening  under  pressure.  Plastic  deformation  in 
TATB  was  investigated  using  Generalized  Stacking  Fault  Energy  (GSFE) 
calculations  at  pressures  of  1  atm  and  5  GPa.  GSFE  calculations,  on  basal  planes, 
predict  that  two  glide  plane  types  exist  for  the  same  glide  plane  normal.  The 
unstable  stacking  fault  energies  on  the  basal  plane  are  less  than  20  mJ/m  at  1  atm, 
indicating  easy  dislocation  glide.  A  stable  stacking  fault  and  a  compound  twin  are 
predicted  to  be  stable  at  1  atm  and  5  GPa. 

Characterization  of  incipient  plasticity  from  nano-indentation  simulations: 

Mechanical  response  of  TATB  during  nano-indentation  was  investigated"  using 
million-atom  MD  simulations.  For  the  (001)  plane,  the  MD  prediction  of  the 
elastic  part  of  the  force-displacement  curve  agrees  well  with  predictions  from 
Hertzian  contact  theory,  obtained  using  the  Stroh  formalism  for  anisotropic 
elasticity  (see,  Fig.  2(a)).  Plastic  deformation  initiates  below  the  indentation 
surface,  progresses  with  kinking  and  delamination  of  basal  planes,  and  results  in 
extensive  pile-up  (see.  Figs.  2(b)  and  (c)),  with  very  high  temperatures  in  the 
vicinity  of  the  indenter  (see.  Fig.  2(d)).  For  the  (100)  plane,  MD  simulations 
predict  an  anomalous,  softer  response  when  compared  to  the  prediction  from 
Hertzian  contact  theory.  (Published  during  second  year  of  project.) 

Determining  the  thermal  conductivity  of  defect  free  TATB  crystals:  We  made 
the  first  predictions  of  the  anisotropic  thermal  conductivity  for  TATB  single 
crystals  using  reverse  non-equilibrium  molecular  dynamics  (rNEMD).'  Our 
conductivity  predictions  were  determined  to  be  generally  insensitive  to  the  details 
of  the  intramolecular  potential,  with  all  variations  predicted  to  be  <12%.^  These 
predictions  are  similar  in  magnitude  to  the  available  experimental  values 
determined  for  pressed  powder  samples.  The  thermal  conductivity  for  TATB  is 
predicted  to  be  considerably  higher  and  more  anisotropic  than  that  for  related 
organic  crystalline  explosives.  Finite-size  effects  were  investigated  and  found  to 
be  minimal,'’^  indicating  that  nanoscale  thermal  transport  in  TATB  crystals  may 
be  well-described  by  continuum  theories  such  as  the  diffusive  heat  equation. 

Understanding  conduction  in  real  (defective)  TATB  crystals:  The  impact  of 
dynamic  structural  transitions  and  molecular  vacancy  defects  on  TATB  crystal 
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thermal  conductivity  was  characterized^  as  a  first  step  toward  understanding  the 
conduction  properties  of  real  (imperfect)  TATB  crystals.  Approximately  linear 
and  anisotropic  relationships  between  molecular  vacancy  defect  density  and 
thermal  conductivity  are  predicted. 

Characterizing  the  (J,  P)  dependence  of  the  conductivity:  TATB  crystal 
thermal  conductivity  was  determined  as  a  function  of  T  and  P  (see  Fig.  3).^  The 
conductivity  is  predicted  to  be  proportional  to  HT,  which  is  consistent  with 
general  analytical  theories,  and  anisotropy  is  predicted  to  decrease  significantly 
with  increasing  T.  Anisotropic  linear  increases  with  increasing  P  are  also 
predicted.  (Published  in  the  second  year.) 

Nanoscale  thermal  transport  in  TATB  is  diffusive:  Molecular  dynamics 
predictions  for  transient  nanoscale  hot  spot  relaxation  were  compared  to 
continuum-based  diffusive  heat  equation  models  and  reveal  such  processes  are 
well-described  by  continuum  approaches  even  at  the  nanometer  length  scale  (see 
Fig.  4).  (Published  in  the  second  year.) 

Anisotropy  in  surface-initiated  melting:  Surface-initiated  melting  was 
investigated”^  using  molecular  dynamics  (MD)  simulations  for  the  three  principal 
crystallographic  planes  exposed  to  vacuum,  with  the  normal  vectors  to  the  planes 
given  by  ftxc,  c^a,  and  axb,  denoted  as  (100),  (010),  and  (001),  respectively.  The 
nature  and  extent  of  disordering  of  the  crystal-vacuum  interface  depends  on  the 
exposed  crystallographic  face,  with  the  (001)  face  exhibiting  incomplete  melting 
and  superheating.  This  is  attributed  to  the  anisotropy  of  the  intermolecular 
hydrogen  bonding  and  the  propensity  of  the  crystal  to  form  stacking  faults  in 
directions  approximately  perpendicular  to  the  (100)  and  (010)  faces.  The  best 
estimate  of  the  melting  temperature  (Tm)  at  1  atm  for  TATB,  obtained  from 
surface-initiated  melting  on  the  (100)  face,  is  851  ±  5  K  and  is  17.7%  higher 
compared  to  the  experimental  value  of  723  K,  obtained  using  pressed  powder 
pellets. 


Calculation  of  the  solid-liquid  coexistence  curve:  Solid-liquid  coexistence 
curve  for  TATB,  at  P  =  1  atm  to  20  kbar,  was  determined^  from  isobaric- 
isothermal  (NPT)  MD  simulations  of  melt  in  contact  with  (100)  face  of  the 
crystal.  An  upper-bound  and  a  lower-bound  for  Tm  at  a  particular  P  is  obtained  by 
performing  these  (solid-liquid  coexistence)  simulations  at  multiple  T  at  a 
particular  P  and  monitoring  the  volume  of  the  simulation  cell  as  a  function  of 
time;  the  system  was  deemed  to  be  above  or  below  Tm  if  the  volume  increased  or 
decreased,  respectively,  with  time.  The  solid-liquid  coexistence  curve  is  obtained 
by  fitting  the  Simon-Glatzel  equation  to  (the  upper-bound)  MD  predictions  (see. 
Fig.  5)  and  is  given  by: 


T  = 

m 


(P+ 6.4808) 
1.6276x10-' 


I  1/ 

72.5892 
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Here,  Tm  is  the  melting  temperature  expressed  in  K  and  P  is  the  pressure 
expressed  in  kbar.  We  note  here  that  the  prediction  of  Tm  at  1  atm,  obtained  from 
coexistence  simulations,  agrees  well  with  the  predictions  from  simulations  of 
surface-initiated  melting. 

Calculation  of  liquid-phase  transport  properties  for  meso-scale  simulations: 

Shear  viscosity  (//),  self-diffusion  coefficient  (Z)),  and  thermal  conductivity  (A)  of 
liquid  TATB,  at  P  =  1  atm  to  20  kbar  and  multiple  temperatures,  are  currently 

O 

being  calculated  from  4  ns  long  isochoric-isothermal  (NVT)  MD  simulations. 
Equilibrium  densities  for  NVT  simulations  are  obtained  from  independent,  2ns 
long,  NPT  simulations.  Table  I  shows  the  values  for  tj,  D,  and  X  obtained  for  Z*  =  1 
atm.  Arrhenius  fit  of  rj  and  D  for  P  =  I  atm  (shown  in  Fig.  6)  predicts  activation 
energy  of  Er,  =  7.29  kcal  mof'  and  Ed  =  12.27  kcal  mof'.  The  thermal 
conductivity  is  predicted  to  be  up  to  80%  lower  than  for  the  crystal  and  linearly 
decreases  with  increasing  T. 

SECOND  YEAR  WORK  (EXCLUDING  ITEMS  MENTIONED  ABOVE 
THAT  CARRIED  OVER  INTO  YEAR  2) 

P-HMX  second-order  elastic  coefficients  on  the  cold  curve:  Molecular 
mechanics  energy  minimization  methods  were  used  to  calculation  to  full  second- 
order  elastic  coefficients  of  P-HMX  as  functions  of  pressure  at  0  K.  From  these, 
the  Voigt-  and  Reuss-average  isotropic  moduli  can  also  be  obtained  as  functions 
of  pressure.  The  results  are  shown  in  Fig.  7. 

Nanoindentation  in  P-HMX:  We  are  using  controlled-displacement  indentation 
methods  similar  to  those  described  above  for  TATB  to  study  the  nanoindentation 
response  for  HMX  surfaces.  As  of  this  writing,  the  results  yield  predictions  in 
close  agreement  with  the  Hertzian  contact  model  predictions  for  low  levels  of 
displacement  but,  in  contrast  to  the  results  for  TATB,  a  lack  of  the  clear  load 
drops  that  would  signify  well-defined  plastic  deformation  events.  Rather,  for  the 
specific  case  of  nanoindentation  normal  to  (110)  we  find  a  gradual  rolloff  behond 
which  the  load-displacement  curve  appears  to  be  flat.  This  led  us  to  investigate  a 
detail  in  our  simulation  protocol,  namely  whether  the  indenter  was  displaced  into 
the  sample  in  a  step-wise  fashion  (using  1.0  A  incremental  displacement  at  2.5  ps 
intervals  during  the  simulation)  versus  in  a  continuous  fusion  (a  very  small 
displacement  at  every  time  step  of  the  simulation  so  as  to  yield  the  same  40  m  s  ' 
overall  indentation  rate  as  used  in  the  stepwise  simulations  and  in  our  published 
paper  for  TATB).  We  find  that  the  general  phenomenon  -  the  smooth  rollover  to  a 
flat  load-displacement  curve  -  is  obtained  in  both  cases  but  (as  one  might  expect) 
the  displacement  at  which  the  rollover  occurs  is  higher  for  the  continuous 
displacement  protocol.  These  results,  which  are  preliminary,  are  shown  in  Fig.  8. 

Kaptiza  resistance:  We  used  the  GCCM  method  to  build  a  (100)/(001)  grain 
boundary  interface  in  TATB  within  an  overall  3-D  periodic  simulation  cell.  The 
grain  boundary  was  oriented  normal  to  the  overall  simulation  cell  and  a  thermal 
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gradient  across  the  grain  boundary  was  prepared  by  defining  “hot”  and  “cold” 
reservoirs  in  the  (001)  and  (100)  regions,  respectively,  far  from  the  grain 
boundary.  The  system  was  simulated  until  it  reached  a  steady-fluctuating  state 
afterwhich  statistics  were  accumulated  for  4  ns.  (See  Fig.  9.)  This  particular  grain 
boundary  orientation  was  chosen  because  it  is  essentially  the  “most  anisotropic” 
one  that  could  exist  between  two  TATB  crystals,  as  is  clear  from  the  figure. 
Despite  the  pronounced  structural  anisotropy,  the  temperature  drop  across  the 
interface  is  only  4  ±  9  K  (based  on  linear  regression  of  temperature  profiles  away 
from  the  reservoirs  or  the  interface),  yielding  a  Kapitza  resistance  of  ~10'  m  K 
W  ',  which  is  orders  of  magnitude  smaller  than  values  obtained  for  covalent 
materials.  Our  initial  explanation  for  this  lack  of  an  effect  for  the  studied  case  is 
that  the  phonon  mean  free  path  length  in  TATB  (estimated  by  us)  and  other 
energetic  crystals  notably  including  RDX  (estimated  by  others)  is  of  the  order  of 
molecular  dimensions  rather  than  the  very  large  values  characteristic  of  covalent 
or  ionic  materials  such  as  diamond,  silicon,  or  UO2. 

TABLES 


Table  I.  Transport  properties  of  liquid  TATB  at  P  =  1  atm. 


T 

(K) 

P 

(kgm-') 

(lO'^Pa  s) 

(10'^m%-') 

2 

(W  m-'  K-') 

900 

1352  (13) 

1.262  (0.006) 

1.0787  (0.0002) 

0.26  (0.01) 

950 

1303  (13) 

1.042  (0.002) 

1.608  (0.001) 

0.25  (0.01) 

1000 

1258 (16) 

0.860  (0.002) 

2.164  (0.001) 

0.23  (0.01) 

1050 

1193  (17) 

0.764  (0.006) 

2.871(0.002) 

0.20  (0.01) 

1100 

1132(19) 

0.640  (0.006) 

3.650  (0.003) 

0.190  (0.009) 

1150 

1046  (23) 

0.497  (0.002) 

5.020  (0.006) 

0.172  (0.008) 
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FIGURES 


Fig.  1. 


TATB  crystal  packing  structure  for  a  3ax3bx4c  supercell  (a)  within  the  planar  molecular 
layers  and  (b)  between  the  molecular  layers.  The  unit  eell  eontains  two  molecules  and 
contributes  two  layers.  Due  to  the  triclinic  symmetry,  the  lattiee  veetors  a,  b,  and  c  are 
not  mutually  orthogonal  so  vectors  a  and  b  in  panel  (a)  and  c  in  panel  (b)  do  not  lie 
exactly  within  the  plane  of  the  page.  Thermal  eonduetion  is  investigated  along  direetions 
a’  and  c'.  Atom  colors  are  gray  for  carbon,  blue  for  nitrogen,  red  for  oxygen,  and  white 
for  hydrogen. 


Nanoindentation  on  the  (001)  plane  of  TATB.  (a)  Foree-displacement  curve  (b) 
Distribution  of  maximum  shear-stress  at  (5  =  9  A  (e)  Pile-up  during  nanoindentation  (d) 
Loeally  averaged  temperature  under  the  nanoindenter. 
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Fig.  3. 


(a)  (b) 


Temperature  (K) 


Anisotropic  thermal  conductivity  X  for  TATB  crystal  (a)  at  P  =  \  atm  and  200  K  <  T  < 
700  K  and  (b)  at  T  =  300  K  and  1  atm  <  P  <25  kbar.  Diamond  and  square  symbols  are 
predictions  from  rNEMD  simulations,  solid  curves  are  weighted  least-squares  regressions 
of  rNEMD  data  fit  to  2  =  mIT  -i-  b  for  (a)  and  to  2  =  mP  +  b  for  (b),  and  dashed  curves  are 
extrapolations  from  the  regressions.  The  inset  in  panel  (a)  shows  the  same  results  as  the 

1  O 

main  figure  but  plotted  against  \IT.  Experimental  values  for  pressed  powder  samples 
are  indicated  by  circle  symbols. 


Fig.  4. 


Temperature  profiles  for  selected  TATB  crystal  layers  in  a  relaxation  simulation  for  an 
idealized  1 -dimensional  hot  spot  that  is  only  three  unit  cells  thick.  Molecular  dynamics 
predictions  (black  curves)  fit  to  a  constant-conductivity  solution  (red  curves)  of  the 
diffusive  heat  equation  (DHE)  are  in  general  qualitative  agreement.  Comparison  to  a 
solution  with  a  temperature-dependent  conductivity  function  (cyan  curves)  highlights  the 
importance  of  including  temperature-dependent  parameters  in  continuum  models. 
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Fig.  5 


Melt-curve  for  TATB  from  solid-liquid  coexistence  simulations. 


Fig.  6. 
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Arrhenius  fits  to  self  -diffusion  coefficient  and  shear  viscosity  of  TATB  at  1  atm. 
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Fig.  7. 
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Second-order  elastic  coefficients  and  derived  isotropic  moduli  for  P-HMX  on  the  cold 
curve.  Only  the  compressive  and  pure-shear  elements  of  the  tensor  are  shown.  Note  the 
small  difference  between  the  Reuss  (uniform  stress)  and  Voigt  (uniform  strain)  results  for 
the  bulk  modulus  in  contrast  to  the  sizeable  difference  for  the  case  of  the  shear  modulus. 
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Fig.  8. 


Nanoindentation  on  (110)  planes  in  P-HMX.  Left:  Simulation  setup.  The  simulation  cell 
is  periodic  in  the  directions  normal  to  the  indentation  directions.  Center:  Load- 
displacement  curves  for  the  step-wise  indentation  (blue  solid  curve)  and  continuous 
indentation  (purple  curve  for  loading  and  green  curve  for  unloading  from  30  A).  The 
Hertzian  prediction  is  shown  as  the  smooth  curve.  Right:  All-atom  snapshot  of  the  system 
for  the  continuous  indenter  displacement  case  for  indentation  depth  30  A. 

Fig.  9. 


Kapitza  resistance  across  a  (001)|(100)  grain  boundary  in  TATB.  The  simulations  are  3-D 
periodic.  A  hot  reservoir  in  the  (001)  region  and  cold  reservoir  in  the  (100)  region  are 
used  to  establish  a  thermal  gradient  in  the  system  that  crosses  the  grain  boundary 
interface.  The  Kapitza  resistance  (thermal  interface  resistance)  is  obtained  by  dividing  the 
temperature  jump  across  the  interface  (obtained  by  linear  extrapolation  of  temperature 
profdes  in  the  respective  phases  far  from  the  reservoirs  and  interface)  by  the  thermal  flux 
(which  is  known  exactly).  Despite  the  extremely  anisotropic  structure  at  the  grain 
boundary,  the  Kapitza  resistance  is  predicted  to  be  only  ~10-11  m  K  W‘  .  This  is 
attributed  to  the  very  small  phonon  mean  free  path  length  characteristic  of  TATB  and 
(apparently)  other  molecular  energetic  crystals. 
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STUDENT  AND  POSTDOCTORAL  RESEARCHER  OUTCOMES 

•  Dr.  Jeffrey  D.  Veals  (MU  postdoctoral  researcher  on  this  project  during  2014)  was 
hired  as  a  postdoc  at  the  Army  Research  Laboratory  -  Aberdeen  Proving  Grounds  in 
September  2014.  He  is  now  an  ARL  staff  member. 

•  Dr.  Matthew  P.  Kroonblawd  (Ph.D.  earned  under  Prof  Sewell,  conferred  May  2016) 
is  now  a  postdoctoral  research  associate  in  the  Physical  and  Life  Sciences  Directorate 
at  Lawrence  Livermore  National  Laboratory. 

•  Dr.  Nithin  Mathew  (MU  postdoctoral  researcher  on  this  project  2013-2016)  has 
accepted  a  postdoctoral  research  associate  position  in  the  Theoretical  Division  at  Los 
Alamos  National  Laboratory  (start  date:  late  2016  or  early  2017). 

•  Dr.  Shan  Jiang  (MU  postdoctoral  researcher  on  this  project  2016)  has  begun  a  tenure- 
track  assistant  professor  position  in  the  Department  of  Mechanical  Engineering  at  the 
University  of  Mississippi  (“Ole  Miss”;  start  date:  November  2016). 

JOURNAL  PUBLICATIONS 

•  Theoretical  determination  of  anisotropic  thermal  conductivity  for  crystalline  1,3,5- 
triamino-2, 4, 6-trinitrobenzene  (TATB),  Matthew  P.  Kroonblawd  and  Thomas  D. 
Sewell,  J.  Chem.  Phys.  139,  074503  (2013). 

•  Theoretical  determination  of  anisotropic  thermal  conductivity  for  initially  defect-free 
and  defective  TATB  single  crystals,  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell, 
J.  Chem.  Phys.  I4I,  184501  (2014). 

•  Generalised  stacking  fault  energies  in  the  basal  plane  of  triclinic  molecular  crystal 

I, 3, 5 -triamino- 2, 4, 6-trinitrobenzene  (TATB),  Nithin  Mathew  and  Thomas  D.  Sewell, 
Philosophical  Magazine  94,  424  (2015). 

•  Strategies  for  non-uniform  sampling  of  molecular  dynamics  phase  space  trajectories 
of  relaxation  phenomena,  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell,  Computer 
Physics  Communications  196,  143  (2015). 

•  Anisotropy  in  surface-initiated  melting  of  the  triclinic  molecular  crystal  1,3,5- 
triamino-2, 4, 6-trinitrobenzene:  A  molecular  dynamics  study,  Nithin  Mathew,  Thomas 
D.  Sewell,  and  Donald  L.  Thompson,  J.  Chem.  Phys.  143,  094706  (2015). 

•  Predicted  anisotropic  thermal  conductivity  for  crystalline  1,3, 5 -triamino-2, 4,6- 
trinitobenzene  (TATB):  Temperature  and  pressure  dependence  and  sensitivity  to 
intramolecular  force  field  terms,  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell, 
Propellants,  Explosives,  Pyrotechnics  41,  502  (2016) 

•  Nanoindentation  of  the  triclinic  molecular  crystal  1,3, 5 -triamino-2, 4,6- 
trinitrobenzene:  A  molecular  dynamics  study,  Nithin  Mathew  and  Thomas  D.  Sewell, 

J.  Phys.  Chem.  C  120,  8266  (2016). 

•  Anisotropic  relaxation  of  idealized  hot  spots  in  crystalline  1 ,3, 5-triamino-2, 4,6- 
trinitrobenzene  (TATB),  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell,  J.  Phys. 
Chem.  C  120,  17214(2016). 

PAPERS  IN  PREPARATION  (meaning  a  complete  manuscript  is  being  finalized): 

•  A  density  functional  theory  study  of  Al  +  CO2  reactions,  J.  D.  Veals,  R.  Chitsazi,  Y. 
Shi,  T.  D.  Sewell,  and  D.  L.  Thompson 
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•  Predicted  melt  curve  and  liquid  state  transport  properties  of  TATB  from  molecular 
dynamics  simulations,  Nithin  Mathew,  Matthew  P.  Kroonblawd,  Thomas  D.  Sewell, 
and  Donald  L.  Thompson 

INVITED  PRESENTATIONS 

•  Matthew  P.  Kroonblawd,  Physical  and  Life  Sciences  Directorate,  Lawrence 
Livermore  National  Laboratory,  January  22,  2016 

•  Nithin  Mathew,  Theoretical  Division,  Los  Alamos  National  Laboratory,  June  21, 

2016 

•  Thomas  S.  Sewell,  Applied  Research  Institute,  University  of  Illinois  at  Urbana- 
Champaign,  July  14,  2014 

•  Thomas  D.  Sewell,  2016  Gordon  Conference  on  Energetic  Materials,  June  7,  2016 

•  Thomas  D.  Sewell,  Los  Alamos  National  Laboratory  MaRIE  Workshop,  Probing 
Dynamic  Processes  in  Soft  Materials  Using  Advanced  Light  Sources,  July  27,  2016 
(post  end-date,  but  covering  much  of  the  work  accomplished  during  project  lifetime) 

•  Thomas  D.  Sewell,  2016  Lawrence  Livermore  National  Laboratory  “Mesoscale 
Modeling  of  Explosives  Initiation  Workshop,”  October  13,  2016  (post  end-date,  but 
covering  much  of  the  work  accomplished  during  project  lifetime) 

•  Thomas  D.  Sewell,  To  be  presented  at  a  Los  Alamos  National  Laboratory  Workshop 
“The  Behavior,  Fabrication  and  Promise  of  Intentionally  Structured  Energetics”,  to  be 
held  10-12  January  2017  (post  end-date,  but  covering  much  of  the  work  accomplished 
during  project  lifetime) 

AWARDS 

•  Matthew  P.  Kroonblawd:  University  of  Missouri-Columbia  Chemistry  Department, 
Breckenridge/Lyons  Award  for  Outstanding  Graduate  Research  (2016) 

•  Thomas  D.  Sewell:  University  of  Missouri-Columbia  Fuldner  Chemistry  Faculty 
Fellow  (2015) 

•  Thomas  D.  Sewell:  University  of  Missouri-Columbia  Fuldner  Chemistry  Faculty 
Fellow  (2016) 

•  Donald  L.  Thompson:  University  of  Missouri-Columbia  Fuldner  Chemistry  Faculty 
Fellow  (2016) 
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University  of  Illinois  at  Urbana-Champaign  (UIUC) 


1.  OVERVIEW  OE  UIUC  SUB-PROJECT 

The  task  objectives  state  in  the  original  proposal  to  be  carried  out  by  Stewart  and  Matalon 
at  Illinois  were  aimed  at  continuum-based  modeling  and  simulation  of  multifunctional 
energetic  materials,  with  the  following  tasks:  1)  Continuum  modeling  and  simulation  of 
microscopic,  multi-phase,  reactive  processes  of  reactants  at  the  nano-  and  micro-scale. 

2)  Design  of  meso-scale,  reactant  micro-structure  to  control  and  fuel-oxidizer  mixing  and 
energy  release-rate  and  3)  System-level  and  performance  evaluation  of  materials. 

A  particular  emphasis  was  placed  on  new  models  relevant  to  condensed  phase 
combustion,  models  that  allow  for  both  solid  and  liquid  and  gas  phase  components,  phase 
change,  mass  diffusion  and  reaction.  The  UIUC  models  are  developed  to  accept  inputs 
from  molecular  dynamic  transport  simulations  carried  out  by  University  of  Missouri  and 
other  researchers,  and  suggest  new  directions  for  simulations. 

At  the  time  of  this  writing  our  efforts  have  completed  2  years  of  funding  and  we  are 
starting  a  3rd  year  of  partial  funding  that  ends  6/14/2015. 

2.  UIUC  ACTIVITIES  DURING  EIRST-YEAR  PERIOD 

A)  Continuum  Model  Development:  Multi-component  and  Mass  Diffusion 

i)  The  Gibbs  formulation:  To  describe  condensed  phase  energetic  material  combustion 
that  occurs  in  both  rocket  propellant  and  liquid  rocket  propellant  ,using  physical  first 
principles,  Illinois  has  come  up  with  a  new  fundamental  continuum  formulation  for  the 
dynamics  of  multi-component,  multi-phase  chemically  reactive  materials.  Since  the 
constituent  materials  used  in  propellants  often  are  used  in  explosives,  and  other  energetic 
materials  such  as  thermite  or  intermetallic  reactive  materials,  we  have  been  able  to  carry 
out  model  development  for  this  AFOSR  effort  that  has  leveraged  off  of  our  work  from 
some  previous  and  ongoing  projects. 

In  particular,  we  note  that  the  Gibbs  formulation  is  a  non-equilibrium,  continuum 
thermodynamic  formulation  that  uses  the  mass  fraction  of  components  (in  all  phases)  as 
the  fundamental  confisuration  variables  instead  of  previously  used  order  parameters  as  a 
marker  of  the  different  phase  field  variables.  Kinetic  equations  that  account  for  the  mass 
creation/destruction  and  transport  of  the  various  components  (species)  control  both  phase 
transformation  (melting  or  evaporation)  and  chemical  reaction.  The  input  to  the  model 
includes  specifying  the  various  components  associated  with  a  given  mixture  of  materials, 
the  equilibrium  (Gibbs)  potentials  for  the  components,  a  single  component  of  stress,  and 
a  single  temperature  in  the  material.  Despite  the  complexity  of  this  model,  it  enables  a 
proper  description  of  many  of  the  effects  one  needs  to  consider  in  the  condensed  phase 
decomposition  and  reaction  processes  of  propellant  materials  at  the  nano  and  micro-scale, 
without  invoking  ad-hoc  parameters.  The  model  allows  for  anisotropic  properties  of 
propellant  crystallites  by  invoking  an  anisotropic  potential  for  the  elastic  solid.  It  also 
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allow  for  anisotropic  thermal  conductivities  and  thermal  expansion,  the  inclusion  of  state 
dependent  reaction  and  diffusion  kinetics  and  the  use  of  Arrhenius  forms  for  chemical 
reaction  rates.  A  significant  amount  of  effort  has  been  expanded  in  order  to  validate  the 
model  and  to  ensure  it  can  solve  the  class  of  propellant  problem  of  interests. 

The  general  Gibbs  formulation  was  completed  this  year  and  is  current  recorded  in  a 
proceedings  paper,  by  Stewart  [5],  and  a  longer  detailed  paper  in  preparation,  that 
includes  many  worked  examples,  see  Stewart  [8]. 

ii)  Multi-component  Diffusion  and  Kinetics  Modeling:  One  aspect  of  our  effort  in  this 
context  has  been  to  include  realistic  multicomponent  diffusion,  as  relevant  to  the 
condensed  material  used  in  propellants,  in  order  to  examine  their  influence  on  the  burning 
(diffusion  flames),  and  energy  released  that  takes  place  at  the  interface  between  the  fuel 
and  oxidizer.  A  series  of  representative  problems  for  condensed  phase  materials  that 
embrace  increasing  complexity  have  been  considered.  The  counterflow  diffusion  flame 
(with  opposed  reactants  meeting  at  a  stagnation  plane)  was  used,  being  a  basic 
configuration  that  represents  the  flow  due  to  deformation,  from  compaction  or  local 
heating  and  thermal  expansion  processes,  in  the  microscale  environment  of  composite 
energetic  materials.  The  multi-component  diffusion  description  uses  a  generalized  Pick 
formulation  with  coefficient  related  to  the  measurable  binary  diffusivities  between  two 
species.  A  fairly  simple  depletion  form  with  Arrhenius  temperature  dependent 
coefficients  is  used  to  describe  the  reaction  rate.  Several  types  of  analyses  were  carried 
out  at  increasing  levels  of  complexities  (see  publications  [1]  and  [3])  and  include:  an 
asymptotic  analysis  valid  in  the  limit  of  low  strain  rates  (high  residence  time  in  the 
reaction  zone),  a  constant  mixture  density  assumption  that  simplifies  the  flow  description, 
diffusion  models  with  equal  and  unequal  molecular  weights  for  the  various  species,  and  a 
full  numerical  study  for  finite  rate  chemistry,  composition-dependent  density  and  strain 
rates  extending  from  low  to  moderate  values.  The  results  suggest  that  in  composite 
materials  there  will  be  a  distribution  of  local  reaction  sites  and  local  frozen  sites  that 
depend  on  the  local  strain  rate.  A  related  (DTRA)  funded  effort  that  focused  at  the  high- 
temperature  titanium/boron  forming  titanium-boride  led  to  a  paper  accepted  for 
publication  in  Combustion  and  Flame  [1]  dealing  with  this  particular  reaction  as  a 
prototype  reaction  (primarily  because  of  the  available  kinetic  data).  Information  collected 
regarding  boron  properties  for  modeling  is  directly  relevant  to  boron  additives  in 
propellants.  Moreover,  the  effort  established  the  basic  methodologies  that  are  needed  to 
analyze  more  complex,  multi-component  kinetics  encountered  in  realistic  propellant 
applications. 

a)  Numerical  verification  studies:  In  addressing  the  aforementioned  studies,  we  noted 
that  the  numerical  problems  that  need  to  be  addressed  are  extremely  stiff,  primarily 
because  of  the  highly  temperature-dependent  diffusivities  and  chemical  reaction  rates, 
and  requires  carrying  out  test  problems  for  verification  purposes.  In  order  to  establish  our 
basic  methods  for  these  very  nonlinear,  multi-scale  problems,  we  carried  out  an  auxiliary, 
verification  study  of  variable  density  saseous  counterflow  flames  131  that  had  not 
previously  been  considered  in  the  literature.  This  work  will  be  submitted  for  presentation 
at  the  next  International  Symposium  on  Combustion  scheduled  in  2016.  We  note  that 
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counterflow  diffusion  flames  are  prototype  flames  used  in  flamelet  modeling  and  relevant 
to  the  non-premixed  environment  associated  with  propellant  applications. 

This  work  was  carried  out  simultaneously  with  that  reported  in  [1]  in  order  to  ensure  that 
our  numerical  tools  and  methodologies  for  these  highly  nonlinear,  multi-scale  problems 
are  sufficiently  robust,  as  we  move  to  model  more  realistic  and  complex  kinetics  required 
for  propellants.  Likewise  in  [4]  Matalon  and  his  student  analyzed  the  structure  and 
dynamics  of  edge  flames  in  shear  flow  and  the  conditions  leading flame  attachment  at 
and  liftoff from  the  base  of  the  injector,  or  flame  blowout.  This  knowledge  is  essential  to 
properly  simulate  diffusion  flames  in  the  non-premixed  environment  of  propellant 
applications. 

b)  Decomposition  schemes  for  molecular  mono-propellants  (explosive).  Crystallites  of 
nitramines  like  HMX  and  RDX  have  been  considered  as  constituents  for  advanced 
propellants  due  to  their  high  energy  density.  These  materials  are  both  monopropellants 
and  explosives.  Their  energy  release  process  involves  both  phase  change,  (melting) 
followed  by  energetic  release  that  occurs  from  reaction  between  molecular  components 
that  evolve  from  the  thermal  decomposition  of  the  crystallites.  Fig.  1  shows  the 
comparison  of  a  Gibbs-based  continuum  simulation  (solid  lines)  of  thermal  explosion  in  a 
constant  volume  process  for  RDX,  compared  with  a  reactive  molecular  dynamics 
simulation  (points)  that  uses  ReaxFF  (carried  out  by  our  UIUC  colleague  S.  Chaudhuri 
and  K  Joshi).  Although  this  work  was  mainly  funded  by  an  complementary  ONR  project, 
it  is  mentioned  here  because  it  is  a  basic  validation  of  the  Gibbs  model  and  is  an  essential 
test  of  the  UIUC  modeling  strategy. 


Mass  fractions 


Pressure 


Temperature 


A  A 

\ 

Figure  1.  Mass  fraction:  liquid  RDX(blue),  intermediate(  red),  gas  phase(  green) 

In  this  case  we  tested  the  concept  of  using  mass  fractions,  defined  by  binning  ranges  of 
molecular  weights  that  are  observed  in  the  reactive  molecular  dynamics  simulation,  and 

assigning  an  average  molecular  weight  for  that  component  in  a  reduced  model.  A  short 

version  of  this  work  will  appear  as  a  reference  conference  proceeding  [8],  and  full  length 
paper  [11],  will  be  submitted  to  the  Journal  of  Chemical  Physics.  This  work  was  the  first 
of  its  kind  and  shows  that  we  can  directly  model  the  results  obtained  from  molecular 

dynamic  simulation  by  continuum  kinetics.  This  was  a  test  of  the  basic  premise  of  our 
Gibbs  formulation. 
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c)  Kinetics  scheme  and  mass  diffusion  for  separated  reactants  in  condensed  phase 

materials:  The  simple  one  step  global,  irreversible  reaetion  F  +  O  — >  P  is  not  suffieiently 
adequate  to  deseribe  realistie  condensed  phase  combustion  processes.  The  burning  of 
aluminum  in  the  presence  of  a  decomposing  oxidizer,  like  Ammonium  Percholorate,  is 
such  an  example.  It  requires  the  consideration  of  more  than  just  three  components  (fuel, 
oxidizer  and  products),  and  include  aluminum  in  its  solid  and  liquid  phases,  the  Oxidizer 
(nominally  AP),  the  decomposition  products  of  AP  (such  as  Perchlorate  liquid),  and  the 
final  products  such  as  aluminum  oxide  in  the  liquid  and  solid  phases.  The  simplest 
scheme  that  retains  the  experimental  facts  suggests  following  the  evolution  of  as  many  as 
8  to  12  components.  If  one  considers  a  one-dimensional  unsteady  problem,  such  as  the 
ignition  at  an  interface,  or  a  two-dimensional  steady  problem  such  as  the  counterflow 
problem,  both  with  reaction  occurring  at  or  near  a  separated  reactants  material  interface, 
then  one  must  provide  a  model  of  mass  diffusion  for  all  the  8  to  12  identified 
components. 

In  January  2015,  at  the  AFRL  Basic  Research  Review  of  Computational  Research,  and 
Nano-Energetics  Review,  held  at  UCLA,  D.  S.  Stewart  presented  a  sketch  of  an 
aluminum/ AP  kinetics  scheme  for  initially  separated  reactants.  Fig.  2  shows  the  outline  of 
a  kinetic  scheme  currently  under  consideration  by  our  group. 


Oxidizer/metal  combustion  (AI/AP)  reaction  diffusion 
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Figure  2.  A  skeletal  kinetics  mechanism  for  AP,  aluminum  combustion  of  separated 
reactants.  Eight  or  more  components  are  required  based  on  work  from  Guirao  and 
Williams' 

A  take-away  point  is  that  one  generally  requires  solid,  liquid  and  gas  components  for 
propellant  combustion  processes  in  the  condensed  phase.  The  kinetic  scheme  proposed 
above  includes  melting  of  the  solid  aluminum  and  solid  oxidizer  (AP),  chemical  reaction, 
as  well  as  the  formation  of  both  solid  and  liquid  aluminum  oxide  products.  The  analysis 
of  this  problem  is  underway  and  will  build  on  early  work  for  thermite  reaction,  [7],  [9] 
that  was  illustrated  for  aluminum/copper  oxide.  The  basic  premise  is  similar  in  that  solid 
reactants  are  initially  separated.  For  the  case  of  aluminum  and  copper  oxide  reactants. 


model  of  ammonium  perchlorate  deflagration  between  20  and  100  atm  ,  C.  Guirao  and  F.  A.  WilliamsAIAA 
Journal,  Vol.  9,  No.  7  (1971),  pp.  1345-1356. 
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both  the  aluminum  and  copper  oxide  melt,  and  the  liquids  react  to  generate  liquid  copper 
and  liquid  and  solid  aluminum  oxide.  And  in  the  same  way  one  must  generate  a  diffusion 
model  that  requires  at  least  an  effective  Fickian  diffusion  matrix.  To  carry  out  such  an 
analysis  molecular  dynamics  simulation  can  be  used  to  great  effect  to  simplify  this 
model,  whereby  one  can  rule  out  the  diffusion  of  certain  components  relative  to  others. 
The  initial  work  of  the  Gibbs  formulation  for  the  simpler  aluminum/copper  oxide 
thermite  was  invaluable  (even  though  it  had  8  components!)  in  developing  schemes  to 
reduce  the  diffusion  matrix  to  a  tractable  reduced  model.  This  point  was  discussed  in  D. 
S.  Stewart  presentation,  "Continuum  Modeling  of  Multiphase  Energetic  Materials  at  the 
Nano  and  Microscale"  at  the  Jan.  2015  UCLA  Workshop. 


iii)  Aluminum  particle  modeling  and  electric  effects: 

It  has  been  hypothesized  that  small  aluminum  particles  as  small  as  nano-sized  particle 
might  provide  enhanced  heat  release  when  exposed  to  oxidizer.  In  particular  a  2010  paper 
by  the  M.  Zacharia  group  at  U.  Maryland  has  suggested  that,  in  the  early  phase  of 
ignition,  at  about  the  temperature  at  which  aluminum  melts,  that  there  is  an  enhanced 
diffusion  mechanism  of  aluminum  cations  due  to  an  induced  electric  charge  separation. 
The  charge  separation  occurs  just  prior  to  ignition.  They  modeled  this  phenomenon  with 
reactive  molecular  dynamics,  using  specifically  ReaxFF. 

a)  Electric  on  Diffusion  Flames:  Matalon  and  his  collaborators  [2]  carried  out  a  precursor 
study  where  they  analyze  ionic  transport  through  the  non-premixed  reaction  zone 
structure.  Two  main  mechanisms  were  identified,  namely  kinetic  and  body  force  effects. 
Kinetic  effects  are  highlighted  by  a  drift  velocity  that  depends  on  the  charge  and  direction 
of  the  electric  field  such  that  positively  and  negatively  charged  species  are  drifted  in 
opposing  directions,  increasing  or  decreasing  the  Fickian  diffusion  velocity.  The  second 
mechanism  is  due  to  the  additional  body  force  highlighted  via  the  “ionic  wind”  that 
affects  the  momentum  balance  and  a  Ohmic  heating  effect  that  could  be  positive  or 
negative  depending  on  the  flow  configuration  and  direction  of  the  electric  field.  Due  to 
the  imposed  symmetry,  the  ionic  wind  effect  in  the  droplet  problem  considered  acts  only 
to  modify  the  radial  pressure  gradient  in  the  combustion  field..  Ohmic  heating  on  the 
other  hand  was  observed  to  act  as  a  heat  source  for  the  droplet  setup,  leading  to  a  slight 
increase  in  flame  temperature,  enhancement  in  the  burning  rate,  and  extension  of  the 
flammability  limits. 

b)  Analysis  of  a  multi-phase  material  model  for  an  Al/AlOx  droplet  In  order  to  use  the 
Gibbs  formulation  to  model  combustion  of  nano-sized  aluminum  particles,  Stewart  and 
student  K.  Lee,  developed  a  model  for  a  spherical  aluminum  particle  with  an  aluminum 
oxide  cap.  Preliminary  work  established  the  stress  and  displacement  distributions  within 
the  coated  aluminum  sphere  and  the  transitions  in  those  states  as  the  temperature  was 
increased.  These  preliminary  calculations  were  reported  at  the  15th  International 


^  On  the  role  of  built-in  electric  fields  on  the  ignition  of  oxide  coated  nano-aluminum:  Ion  mobility  versus  Fickian 
diffusion,  Brian  J.  Henz,  Takumi  Hawa,  and  Michael  R.  Zachariah,  J.  Appl.  Phys.  107,  024901  (2010) 
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Numerical  Combustion  meeting,  Spring  2015,  [14], 


This  work  is  the  basis  for  our  undergoing  analysis  of  oxidation  (burning)  of  nano,  micro¬ 
sized  aluminum  (or  metal)  particle  that  occur  by  diffusional  processes.  In  which  case  will 
will  add  an  aluminum  cation  and  and  oxygen  anions  that  will  diffusive  through  the  oxide 
layer,  and  use  an  approach  similar  to  that  Matalon  used  in  [2],  We  expect  to  model 
Zacharia-group  configuration  for  the  particle  with  the  Gibbs  formulation,  through  the  end 
of  this  grant.  And  hope  to  confirm  our  model  results  with  MD  simulations  carried  out  the 
the  Sewell-group  at  U  Missouri  or  the  Chaudhuri  group  here  at  UIUC. 


Publications  and  Presentations  of  the  Stewart/Matalon-Group  2014-2015 

*  -  cites  AFOSR/RTE  support 

Peer  Reviewed  Publications  (accepted  or  submitted) 

[1.]  *  Diffusion  flames  in  condensed-phase  energetic  materials:  Application  to 
Titanium-Boron  combustion,  Sushilkumar  P.  Koundinyan,  John  B.  Bdzil,  Moshe 
Matalon,  and  D.  Scott  Stewart,  to  appear  Combustion  and  Flame  (2015) 

[2.]  *  Electric  Field  Effects  in  the  Presence  of  Chemi-Ionization  on  Droplet  Burning 

A.  Patyala  ,  D.  Kyritsis  ,  M.  Matalon  (submitted  to  Combustion  and  Flame) 

[3.]  *  Symmetric  counterflow  diffusion  flames,  S.  Koundinyan,  M.  Matalon,  D.  S. 
Stewart,  to  be  submitted  to  the  36th  International  Symposium 

[4.]  *  Structure  and  dynamics  of  edge  flames  in  the  near  wake  of  unequal  merging 
shear  flows,  Z.  Lu  and  M.  Matalon,  submitted  to  Combustion  Theory  and  Modelling. 

Peer  Reviewed  Proceedings  Papers 

[5.]  *  Modeling  reaction  fronts  of  separated  condensed  phase  reactants,  Sushil 
Koundinyan,  M.  Matalon,  D.  S.  Stewart,  and  J.  Bdzil,  will  appear  in  Shock  Compression 
of  Condensed  Matter  -  2015:  Proceedings  of  the  American  Physical  Society  Topical 
Group  on  Shock  Compression  of  Condensed  Matter;  AIP  Conf  Proc.,  Tampa,  Florida, 
June  2015. 

[6.]  *  A  Gibbs  Formulation  for  Reactive  Materials  with  Phase  Change  ,  D.  Scott 
Stewart,  will  appear  in  Shock  Compression  of  Condensed  Matter  -  2015:  Proceedings  of 
the  American  Physical  Society  Topical  Group  on  Shock  Compression  of  Condensed 
Matter;  AIP  Conf  Proc.,  Tampa,  Florida,  June  2015. 

[7.]  *  Modeling  the  Shock  Ignition  of  a  Copper  Oxide  Aluminum  Thermite  ,  Kibaek 
Lee,  D.  S.  Stewart,  M.  Clemenson,  N.  Glumac,  C.  Murzyn,  will  appear  in  Shock 
Compression  of  Condensed  Matter  -  2015:  Proceedings  of  the  American  Physical  Society 
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Topical  Group  on  Shock  Compression  of  Condensed  Matter;  AIP  Conf  Proc.,  Tampa, 
Florida,  June  2015. 

[8.]  Mirrored  continuum  and  molecular  scale  simulations  of  the  ignition  of  gamma 
phase  RDX,  D.  Scott  Stewart,  Santanu  Chaudhuri,  Kaushik  Joshi.  Kiabek  Lee,  will 
appear  in  Shock  Compression  of  Condensed  Matter  -  2015:  Proceedings  of  the  American 
Physical  Society  Topical  Group  on  Shock  Compression  of  Condensed  Matter;  AIP  Conf 
Proc.,  Tampa,  Florida,  June  2015. 


Papers  in  Preparation  for  Peer  Reviewed  Publication 
(with  partial  or  complete  typescript) 


[9.]  *  Modeling  the  Shock  Ignition  of  a  Copper  Oxide  Aluminum  Thermite,  Kibaek 
Lee,  D.  S.  Stewart,  to  be  submitted  to  Journal  of  Applied  Physics 

[10.]  *  A  Gibbs  Formulation  for  Reactive  Materials  with  Phase  Change,  D.  Scott 
Stewart,  in  preparation. 

[11.]  Mirrored  continuum  and  molecular  scale  simulations  of  the  ignition  of 
gamma  phase  RDX,  D.  Scott  Stewart,  Santanu  Chaudhuri,  Kaushik  Joshi.  Kiabek  Lee, 
to  be  submitted  to  the  Journal  of  Chemical  Physics. 

Conference  Presentations  without  Proceedings 

[12]  *  Compressible  mixing  of  liquid  fuels,  Alberto  Hernandez,  Svjetlana  Stekovic,  D. 
Scott  Stewart  and  Alexander  Wass,  15th  International  Conference  on  Numerical 
Combustion,  Avignon,  France,  April  2015 

[13]  *  Predicting  extinction  in  condensed  phase  combustion  in  a  counterflow 
geometry,  Sushil  Koundinyan,  M.  Matalon,  D.  S.  Stewart,  and  J.  Bdzil,  15th 
International  Conference  on  Numerical  Combustion,  Avignon,  France,  April  2015 

[14]  *  Analysis  of  a  multi-phase  material  model  for  an  Al/AlOx  droplet,  Kibaek  Lee, 
D.  Scott  Stewart,  15th  International  Conference  on  Numerical  Combustion,  Avignon, 
France,  April  2015 
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3. 


UIUC  ACTIVITIES  DURING  SECOND- YEAR  PERIOD 


A)  Continuum  Model  Development:  The  Gibbs  formulation:  Modeling  of  HMX 
decomposition. 

To  describe  condensed  phase  energetic  material  combustion  that  occurs  in  both  rocket 
propellant  and  liquid  rocket  propellant,  using  physical  first  principles,  Illinois  developed 
a  new  fundamental  continuum  (Gibbs)  formulation  for  the  dynamics  of  multi-component, 
multi-phase  chemically  reactive  materials.  Since  the  constituent  materials  used  in 
propellants  often  are  used  in  explosives,  and  other  energetic  materials  such  as  thermite  or 
intermetallic  reactive  materials,  we  have  been  able  to  carry  out  model  development  for 
this  AFOSR,  leveraged  off  of  our  work  from  some  previous  and  ongoing  projects.  In 
particular  a  work  recently  published  in  [1],  used  a  Gibbs  formulation  that  and  defined  a 
kinetics  scheme  for  nitramine-based  energetic  materials,  that  was  used  modeled 
decomposition  in  RDX.  The  Gibbs  formulation  faithfully  modeled  the  averages  observed 
in  reactive  molecular  dynamics  simulations  and  thus  defined  relevant  intermediate 
species  and  phases. 

Modeling  HMX  Decomposition  for  Propellants 
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Figure  I.  A  simplified  HMX  decomposition  scheme  for  propellants  that  includes  changes 
from  beta  to  delta  solid  phase,  liquid  phase,  gas  and  product  gas  phase  components. 

Notably  HMX  is  a  constituent  used  in  high  performance  ballistic  missile  solid  rocket 
motors.  A  similar  kinetics  model  was  employed,  and  significantly  extended  to  model  a 
HMX  decomposition  flame.  Figure  1.  displays  the  simplified  kinetics  and  identifies  the 
component  used  in  the  model.  The  UIUC  model  now  includes  transport  and  allows  for 
the  progression  of  a  decomposition  path,  expected  to  be  observed  in  an  HMX 
monopropellant  flame  that  goes  from  solid  phase  through  liquid  to  vapor  to  burnt 
product. 
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Figure  2.  Schematic  of  the  Strand  burner  configuration  for  monopropellant  burning  of 
HMX,  developed  in  the  UM/CIT/UIUC  joint  effort.  The  UIUC  model  developed  last  year 
considered  the  phase  shown  as  basic  decomposition  steps. 

A  ID  unsteady  model  for  burning  of  HMX  with  dissolved  liquids  and  gases. 

Near  the  end  of  the  work  period,  UIUC  had  developed  a  one-dimensional,  but  unsteady 
model,  for  the  burning  of  HMX.  The  model  had  the  multiple  components  discussed 
above  and  allows  for  dissolved  liquids  and  gases  in  the  material  mixture.  The  reduced 
model  allowed  one  to  solve  for  the  thermal  and  mass  fractions  and  then  compute  the 
stress/strain  fields.  The  modeling  framework  is  tractable  and  will  be  passed  to  the  CIT 
(Ortiz)-group  if  there  is  additional  AFOSR  funding  for  this  work.  The  detailed 
derivations  of  the  model  are,  as  of  yet,  unpublished,  but  they  are  extensive  and  the  model 
is  complete  if  not  yet  validated,  [2].  Specifically,  the  intent  is  to  continue  study  the 
properties  of  this  model  and  generate  a  still  reduced  (but  tractable)  model  that  can  be 
implemented  into  the  Ortiz  computational  modeling  framework. 

A  series  of  ID,  test  simulations  were  carried  out  using  the  model.  The  left  edge  of  a  100- 
micron  slab  was  raised  to  constant  temperature  at  the  left  edge  and  insulated  on  the  right 
edge  See  Fig.  3.  Then  the  simulations  showed  phase  transition,  melting  and  incipient 
reaction.  Therefore  we  were  able  to  verify  the  basic  integrity  of  the  model  and  its 
feasibility.  Importantly  the  thermal  expansion  in  the  material  is  not  negligible  and 
significant  stress  fields  can  be  developed  if  the  material  slabs  are  confined  during  thermal 
decomposition.  While  this  is  not  unexpected  we  demonstrated  that  the  model  could 
compute/model  these  affects  that  are  almost  never  accounted  for  in  other  models.  Figure 
4  shows  typical  results  that  illustrate  that  for  a  confined  slab  up  to  2  GPa,  (20  Khar)  of 
pressure  maybe  realized  in  a  confined  slab  of  HMX  propellant  on  the  millisecond  time 
scale. 

At  the  time  of  this  writing  this  work  is  still  in  development  and  will  result  in  publication 
in  the  future.  This  work  can  be  the  basis  for  follow  on  research  supported  by  AFOSR,  and 
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will  be  directed  to  integrated  simulations  with  Cal-Tech. 
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1200  K 


Vcipcjrto 
gas  products 


704  K 

552  K 
438  K 


iki'jidto  v-ipor 

delta  to  liquid 

'  to  delta 


rlj'jJyl 

iil-jb 


insulated 


(Til  =  0 
stress  normal 
to  slab 


rigid  wall  zero 
displacement 


■ .  298.15  K 

I - 

0  100  M 


Figure  3.  Basic  set  up  for  the  test  problem.  The  left  part  of  the  figure  shows  the 
temperatures  and  the  associated  phase  transitions.  The  figure  is  labeled  with  the  boundary 
conditions  imposed  on  a  100  micron  slab  with  a  left  edge  raised  and  kept  at  1200  Kelvin. 
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other  Publications  and  Works  Supported  by  this  Grant. 


During  the  last  year  of  the  grant,  a  number  of  other  works  were  eompleted  and  published 
by  Matalon  and  Stewart  and  their  students  at  UIUC.  The  publieations  that  were 
eompleted  are  listed  below  with  attaehed  abstraets  that  provide  a  summary  of  the  work. 
In  all  eases  the  work  earried  out  was  supported  wholly  or  in  part  by  this  grant  and  eite 
AFOSR  support.  These  work  treat  a  variety  of  important  issued  that  inelude  multi- 
eomponents  modeling,  thermal  transport  and  eleetrie  field  variation  and  their  impaets  on 
flames.  All  these  issues  are  relevant  to  propellant  flame  modeling. 


Diffusion  flames  in  condensed-phase  energetic  materials:  Application  to  Titanium- 
Boron  combustion 

Combustion  and  Flame  Vol.  162  (2015)  pp.  4486-4496 

Sushilkumar  P.  Koundinyan,  John  B.  Bdzil,  Moshe  Matalon  and  D.  Scott  Stewart 

The  characteristics  of  a  steady  diffusion  flame  that  arises  at  the  interfaces  of  two  condensed  phase  reactant 
streams  that  form  an  opposed  counterflow  are  discussed.  We  assume  that  the  flow  is  due  to  deformation 
from  compaction  or  local  heating  and  thermal  expansion  processes  in  the  microscale  environment  of 
composite  energetic  materials.  As  a  representative  example  of  high  temperature  combustion  of 
metal/intermetallic  reactants,  the  overall  reaction  of  titanium  and  boron  to  create  titanium  diboride  products 
is  considered  under  near  isobaric  conditions.  The  multi-component  diffusion  description  uses  a  generalized 
Tick  formulation  with  coefficients  related  to  the  binary  diffusivities  defined  in  the  Maxwell-Stefan 
relations.  A  fairly  simple  depletion  form  with  Arrhenius  temperature  dependent  coefficients  is  used  to 
describe  the  reaction  rate.  Several  types  of  analyses  are  carried  out  at  increasing  levels  of  complexities:  an 
asymptotic  analysis  valid  in  the  limit  of  low  strain  rates  (high  residence  time  in  the  reaction  zone),  a 
constant  mixture  density  assumption  that  simplifies  the  flow  description,  diffusion  models  with  equal  and 
unequal  molecular  weights  for  the  various  species,  and  a  full  numerical  study  for  finite  rate  chemistry, 
composition-dependent  density  and  strain  rates  extending  from  low  to  moderate  values.  All  are  found  to 
agree  remarkably  well  in  describing  the  flame  structure,  the  flame  temperature  and  the  degree  of 
incomplete  combustion.  Of  particular  importance  is  the  determination  of  a  critical  strain  rate  beyond  which 
steady  burning  may  no  longer  be  observed.  The  analysis  has  a  general  character  and  can  be  applied  to  other 
condensed  phase  energetic  material  systems,  where  reaction  and  diffusion  occur  in  the  presence  of  flow  and 
material  deformation. 


Structure  and  dynamics  of  edge  flames  in  the  near  wake  of  unequal  merging  shear 
flows 

Combustion  Theory  and  Modelling,  Vol.  20,  No.  2  (2016)  pp.  258-295 
Zhanbin  Lu  and  Moshe  Matalon 

We  examine  in  this  study  the  structure  and  dynamic  properties  of  an  edge  flame  formed  in  the  near-wake  of 
two  initially  separated  shear  flows,  one  containing  fuel  and  the  other  oxidiser.  A  comprehensive  study  is 
carried  out  within  the  diffusive -thermal  framework  where  the  flow  field,  computed  a-priori,  is  used  for  the 
determination  of  the  combustion  field.  Our  focus  is  on  the  effects  of  three  controlling  parameters:  the 
Damkohler  number  controlling  the  overall  flow  rate,  the  oxidiser-to-fuel  strain  rate  ratio  of  the  supply 
streams  that  determines  the  extent  of  oxidiser  entrainment  towards  the  mixing  zone,  and  the  Lewis  number, 
assumed  equal  for  the  fuel  and  oxidiser,  that  depends  on  the  mixture  composition.  Response  curves, 
representing  the  edge  flame  standoff  distance  as  a  function  of  Damkohler  number,  exhibit  two  distinct 
shapes:  C-shaped  and  U-shaped  curves  characterising  the  response  of  low  and  high  Lewis  number  flames, 
respectively.  Stability  considerations  show  that  the  upper  solution  branch  of  the  C-shaped  response  curve  is 
unstable  and  hence  corresponds  to  physically  unrealistic  states,  but  due  to  heat  conduction  toward  the  cold 
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plate  the  lower  solution  braneh  is  always  stable.  The  states  forming  this  solution  braneh  eorrespond  to 
flame  attaehment,  where  the  edge  flame  remains  praetieally  attaehed  to  the  tip  of  the  plate  until  it  is  blown 
off  by  the  flow  when  the  veloeity  exeeeds  a  eritieal  value.  The  U-shaped  response,  on  the  other  hand, 
eonsists  of  equilibrium  states  that  are  globally  stable.  Thus,  high  Lewis  number  flames  ean  be  always 
stabilised  near  the  splitter  plate,  with  the  edge  held  stationary  or  undergoing  a  baek  and  forth  motion,  or 
lifted  and  stabilised  downstream  by  the  flow.  Insight  into  the  distinet  stabilisation  eharaeteristies,  exhibited 
by  the  different  Lewis  number  eases,  is  given  by  examining  the  relationship  between  the  loeal  flow  veloeity 
and  the  edge  propagation  speed. 

Electric  field  effects  in  the  presence  of  chemi-ionization  on  droplet  burning 

Combustion  and  Flame  Vol.  164  (2016)  pp.  99-110 
Advitya  Patyal,  Dimitrios  Kyritsis  and  Moshe  Matalon 

The  effeets  of  an  externally  applied  eleetrie  field  on  the  burning  eharaeteristies  of  a  spherieally  symmetrie 
fuel  drop,  ineluding  the  flame  strueture,  flame  standoff  distanee,  mass  burning  rate  and  flame  extinetion 
eharaeteristies  of  the  diffusion  flame  are  studied.  A  redueed  three-step  ehemieal  kinetie  meehanism  that 
refleets  the  ehemi-ionization  proeess  for  general  hydroearbon  fuels  has  been  proposed  to  eapture  the 
produetion  and  destruetion  of  ions  inside  the  flame  zone.  Due  to  the  imposed  symmetry,  the  effeet  of  the 
ionie  wind  is  simply  to  modify  the  pressure  field.  Our  study  thus  foeuses  exelusively  on  the  effeets  of 
Ohmie  heating  and  kinetie  effeets  on  the  burning  proeess.  Two  distinguished  limits  of  weak  and  strong 
field  are  identified,  highlighting  the  relative  strength  of  the  internal  eharge  barrier  eompared  to  the 
externally  applied  field.  For  both  limits,  signifieantly  different  eharged  speeies  distributions  are  observed. 
An  inerease  in  the  mass-burning  rate  is  notieed  with  inereasing  the  strength  of  the  eleetrie  field  in  both 
limits,  with  a  small  ehange  in  flame  temperature.  Inereasing  external  voltages  pushes  the  flame  away  from 
the  droplet  and  eauses  a  strengthening  of  the  flame  with  a  reduetion  in  the  extinetion  Damkohler  number. 

Temperature-dependent  Transport  and  Thermal-diffusion  Effects  on  Diffusion 
Flames 

Submitted  to  Combustion  Theory  and  Modelling 
Joseph  E.  Hibdon  Jr.  and  Moshe  Matalon 

The  effeets  of  temperature-dependent  transport  and  thermal  diffusion,  also  known  as  Soret-speeies 
transport,  on  the  strueture  and  eharaeteristies  of  a  diffusion  flame  are  studied.  The  eonfiguration  adopted  is 
a  planar  unstrained  diffusion  flame,  where  a  uniform  flow  eontaining  either  fuel  or  oxidizer  is  direeted 
toward  the  reaetion  zone  with  the  oxidizer  or  fuel  diffusing  from  the  opposite  boundary  against  the  stream. 
Ineluded  in  this  diseussion  is  the  no-flow  ease,  where  the  reaetants  reaeh  the  reaetion  zone  purely  by 
diffusion.  The  model  allows  for  non-unity  and  distinet  Lewis  numbers  for  the  fuel  and  oxidizer  and  for 
finite-ehemistry  effeets,  eovering  the  entire  range  of  Damkohler  numbers,  from  eomplete  eombustion  down 
to  extinetion.  The  asymptotie  analysis  follows  the  earlier  studies  on  diffusion  flames,  extended 
appropriately  to  aeeommodate  for  the  temperature-dependent  transport  eoeffieients  and  diffusion  fluxes 
that  result  from  the  Soret  effeet.  The  results  show  that  both  effeets  lead  to  a  shift  in  the  position  of  the 
reaetion  sheet  with  a  eorresponding  modifieation  of  the  stoiehiometrie  flame  temperature.  Aeeounting  for 
temperature-dependent  transport  eoeffieients  lowers  the  eritieal  Damkohler  number,  below  whieh  flame 
extinetion  oeeurs,  thus  extending  the  flammability  limits  eompared  to  predietions  obtained  when  these 
eoeffieients  are  retained  eonstant.  Soret  effeets  have  also  an  effeet  on  these  flame  eharaeteristies,  but  the 
trend  depends  on  the  fuel  type.  The  standoff  distanee  of  the  reaetion  sheet  shifts  towards  the  fuel/oxidizer 
for  “heavy/light  fuels”,  eausing  a  reduetion/inerease  in  flame  temperature,  respeetively.  The  effeet  on  flame 
extinetion  beeomes  signifieant  only  in  “heavy  fuels”  leading  to  a  reduetion  in  the  extinetion  Damkohler 
number. 

Counterflow  Diffusion  Flames  -  Numerical  Simulations  vs  Asymptotic 

Submitted  to  Combustion  theory  and  Modelling 

Sushilkumar  P.  Koundinyan,  John  B.  Bdzil,  Moshe  Matalon  and  D.  Seott  Stewart 
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Numerical  simulations  of  counterflow  diffusion  flames  that  properly  accounts  for  thermal  expansion  are 
presented  along  with  results  of  the  corresponding  asymptotic  theory  for  large  activation  energies.  Unlike 
the  idealized  constant  density  case,  the  variable  density  asymptotic  formulation  requires  integrating  the 
outer  convective-diffusive  equations  numerically  while  satisfying  the  jump  relations  that  mimic  the 
reaction  diffusion  processes  inside  the  reaction  sheet.  The  numerical  implementation  of  the  asymptotic 
theory,  and  the  numerical  simulations  of  the  full  governing  equations  lead  to  qualitatively  identical  results, 
and  show  good  quantitative  agreement  as  well.  In  particular,  the  extinction  Damkohler  number  and  the 
flame  temperature  and  position  at  extinction  obtained  from  both  formulations  are  in  good  agreement  for  a 
wide  array  of  parameters. 
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California  Institute  of  Technology 


CALTECH  ACTIVITIES  DURING  FIRST-YEAR  PERIOD 

The  work  tasked  to  Caltech  is  concerned  with  the  development  of  mesoscale  simulations  and 
theory  for  multi-material  and  polycrystalline  and  microstructured  materials.  Specific  subtasks  are  the 
generation  of  useable,  verified  and  validated  material  descriptions,  mixtures  and  interfacial  material 
interactions  that  can  be  directly  incorporated  into  any  larger  continuum-based  computational 
framework.  A  primary  focus  of  the  work  specifically  concerns  the  development  of  hierarchical 
theoretical  methods  for  understanding  and  predicting  anisotropic  thermal  transport  and  energy 
release  in  rocket  propellant  formulations.  An  even  more  specific  focus  concerns  the  development  of 
understanding  and  predictive  capability  regarding  anisotropic  thermal  transport  and  energy  release  in 
advanced  rocket  propellants,  leading  to  a  practical  capability  for  a  priori  propellant  design. 
Propellants  have  traditionally  been  modeled  using  hydrocodes  that  neglect  entirely  the  strength  of 
the  material.  However,  recent  work  (Kroonblawd  &  Sewell,  J.  Chem.  Phys.,  2013)  has  shown  that 
the  thermal  conductivity  of  energetic  materials  such  as  TATB  can  be  extremely  anisotropic,  with 
conductivity  contrasts  of  orders  of  magnitude  depending  on  transport  direction.  It  seems  likely  that 
this  extreme  anisotropy  should  in  turn  result  in  a  strong  directional  dependence  of  reaction-front 
speeds,  a  phenomenon  that  appears  to  be  poorly  understood  at  present.  This  acute  sensitivity  to  local 
anisotropy  in  turn  raises  opportunities  for  the  optimal  design  of  microstructures  resulting  in 
enhanced  engineering  properties  of  propellants,  such  as  specific  impulse,  a  design  tradeoff  that  also 
appears  to  be  unexplored  at  present. 

Work  at  Caltech  to  date  has  progressed  on  two  different  fronts,  described  next  in  turn. 

Fundamental  understandins  of  the  effect  of  anisotropy  on  reaction-front  speeds.  The  essential 
difficulty  in  developing  such  understanding  is  the  multiscale  nature  of  the  phenomenon,  which 
includes  both  atomic  level  rate-limiting  processes,  such  as  thermal  vibrations  and  bond  breaking, 
collisions,  coupled  to  slow  processes  such  as  heat  and  mass  transfer,  viscosity  and  complex  reactions 
paths.  Indeed,  the  challenge  of  performing  atomistic  full-chemistry  simulations  over  time-scales 
relevant  to  macroscopic  chemical  processes  and  diffusion-mediated  phenomena  constitutes  a  long¬ 
standing — yet  stubbornly  elusive — goal  in  computational  science.  For  propellants,  the  time-scale  gap 
is  staggering,  from  thermal  vibrations  (femptosecond)  to  device  operation  (seconds).  The  spatial- 
scale  gap  is  equally  staggering,  from  atomic  lattice  scale  (Angstroms)  to  device  dimensions  (m).  And 
yet,  neither  of  the  limiting  scales  is  sufficient  to  provide  a  good-enough  handle  for  material  and 
device  design:  On  one  had  we  wish  atomistic  realism  across  within  reaction  zone  in  order  to  account 
for  complex  difficult-to-predict  reaction  paths;  on  the  other  hand  we  wish  to  make  contact  with 
engineering  systems  and  data.  The  fundamental  question  is,  therefore:  How  to  effect  space- time 
coarse-graining  (atoms  to  device)  without  introducing  spurious  physics  and  without  essential  loss  of 
information?  In  effect,  we  wish  for  atomistic  realism  without  the  need  to  resolve  the  length  scale  of 
the  crystal  lattice  and  the  time  scale  of  thermal  vibrations. 

We  have  addressed  under  the  auspices  of  the  grant  this  challenge  by  developing  a  theory  of 
mesodynamics.  The  theory  is  based  on  a  maximum-entropy  (max-ent)  formulation  of  non¬ 
equilibrium  statistical  mechanics  that  allows  temperature  to  be  defined  locally,  atom  by  atom.  In  this 
manner,  we  achieve  a  statistical  description  of  thermal  vibrations  away  from  equilibrium.  A  system 
of  mesodynamical  equations  is  then  defined  by  projecting  the  Liouville  equation  onto  the  space  of 
max-ent  measures.  The  resulting  equations  reduce  to  Hamilton’s  equations  at  zero  temperature,  in 
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which  limit  the  energy  equation  beeomes 
redundant.  Mesodynamies  essentially 
combines  two  different  deseriptions  of  the 
motion:  deterministic,  for  the  slow 

components  of  the  motion  that  are  not 
resolved  by  the  time  step;  and  statistical,  for 
the  remaining  fast  components  of  the  motion, 
such  as  thermal  vibrations,  not  resolved  by 
the  time  step.  In  caleulations,  we  use 
asynchronous  integrators  that  allow  each 
atom  to  be  updated  on  its  own  time  step,  thus 
providing  various  degrees  of  time  resolution  depending  on  local  resolution  requirements.  In  addition, 
the  meehanical  energy  of  the  fast  eomponents  of  the  motion  not  resolved  by  the  time  step  is 
converted  locally  to  heat  by  recourse  to  stiffness-proportional  damping.  For  harmonie  systems,  this 
form  of  dissipation  damps  out  preferentially  the  high-frequeney  components  of  the  motion  while 
leaving  low-frequency  components  comparatively  unchanged.  Conversely,  when  rare  events  oeeur, 
such  as  hops  or  collisions,  the  local  time  step  is  decreased  in  order  to  provide  fine  temporal 
resolution  of  the  events.  When  such  is  done,  heat  is  pumped  back  into  the  mechanical  equations 
using  a  random-phase  approximation.  In  this  manner,  we  aehieve  a  fiilly-resolved,  deterministic, 
description  of  energetic  rare  events,  such  as  bond-breaking  interaetions  during  ehemieal  reactions, 
while  simultaneously  achieving  a  coarse-grained  statistieal/  mesodynamical  deseription  elsewhere. 
We  have  already  aehieved  accelerations  of  orders-of-magnitude  in  test  eases.  At  present,  we  are 
testing  the  approach  in  applications  concerned  with  the  oxidation  of  magnesium  and  graphite  (figure, 
top  left).  In  addition  to  providing  test  cases  for  the  development  of  the  mesodynamies  solver,  these 
applications  are  expected  to  shed  light  on  the  role  of  thermal  anisotropy  on  reaetion-front  kineties. 


Development  of  methods  to  optimize  propellant  microstructures.  The  aim  of  this  effort  build  on 
fundamental  information  from  atomic-scale  ealeulations,  ineluding  orientation-dependent  reaction- 
front  speeds,  in  order  to  optimize  parameterized  propellant  mierostructures  of  interest,  including 
composites  eomprising  a  Polymer  binder,  oxidizer  granules  and  additives  (e.g.,  NH4C104 
Composite  Propellant,  APCP).  Key  design  parameters  to  be  optimized  include  volume  fractions, 
textures  and  grain/inelusion  morphologies.  We  connect  micromechanical  properties  to  overall 
propellant  performanee  through  front  tracking  simulations  of  the  passage  of  reaction  fronts  through 
representative  volume  elements,  thus  determining  effeetive/maeroscopic  burning  rates.  Such 
effective  burning  rates  are  subsequently  optimized  for  maximum  speeific  impulse.  To  date,  we  have 
completed  the  development  and  testing  of  a  ray-tracing  solver  with  capability  for  tracking  the 
passage  of  reaetion  fronts  through  arbitrary  locally-anisotropic  spatial  mierostruetures  (figure,  next 
page).  At  present,  the  loeal  reaetion-front  speed  depends  on  the  local  anisotropic  thermal 
conductivities  through  a  simple  one-step  reaction  model,  but  the  eomputational  capability  is 
sufficiently  general  to  integrate  micromechanieal  models  as  they  become  available.  The  front 
traeking  ealeulations  are  fast  enough  as  to  make  microstructural  optimization  feasible.  In  addition  to 
these  eomputational  developments,  at  present  we  are  attempting  to  understand  the  multiseale 
structure  of  the  problem  by  analytical  means,  including  closed-form  charaeterization  of  effeetive 
behavior  in  the  sharp-interfaee  limit  and  in  the  homogenization  limit. 


CALTECH  ACTIVITIES  DURING  SECOND-YEAR  PERIOD 
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1  Introduction 

The  specific  Impulse  delivered  by  present-day  propellants  has  remained  flat 
since  the  1 980's.  It  Is  clear  that  radically  new  propellant  preparation  paradigms 
are  needed  In  order  to  break  the  current  impasse.  Recent  advances  In  manu¬ 
facturing,  3D  printing,  advanced  lithography,  nanocomposite  materials,  vastly 
expand  the  design  space  of  propellants  and  call  for  a  systematic  reevaluation 
of  solid  propellant  technology.  Already,  adding  spheroidal  aluminum  particles 
and  other  metallic  phases  to  propellants  Is  a  common  technique  to  increase 
thrust.  However,  such  approaches  beg  the  question  of  optimally  engineered 
propellants.  A  fundamental  question  In  optimal  design  concerns  the  geome¬ 
try  and  topology  of  metallic  phases  resulting  In  the  best  possible  propellant 
performance. 

We  have  developed  capability  that  enables  topology  optimization  of  com¬ 
plex  propellant  designs  that  best  balance  multiple  competing  objectives:  weight, 
durability,  conductivity,  manufacturability  and  others.  We  have  specifically 
applied  the  capability  to  pilot  studies  of  aluminized  composites.  We  assume 
that  the  reaction  rate  of  a  propellant  Is  limited  by  Its  effective  thermal  conduc¬ 
tivity,  as  the  propellant  needs  to  transport  heat  to  the  reaction  front  In  order 
to  sustain  the  reaction.  This  ansatz  suggests  maximizing  effective  thermal 
conductivity  as  one  of  the  optimization  objectives.  However,  thermal  conduc¬ 
tivity  needs  to  be  carefully  balanced  against  competing  objectives  such  as 
structural  Integrity  and  chemistry.  In  what  follow,  the  multi-objective  topol¬ 
ogy  optimization  developed  to  date  is  showcased  by  means  of  test  examples. 
We  anticipate  that.  In  subsequent  phases  of  the  study,  this  capability  will 
supply  an  effective  tool  for  uncovering  novel  engineered  propellant  designs 
with  exceptional  performance. 

2  Review  of  previous  work 

Topology  optimization  with  the  objective  of  maximizing  heat  conduction  en¬ 
ables  the  design  of  structures  that  may  effectively  dissipate  or  transmit  heat 
generated  by  a  source.  This  design  problem  has  been  explored  In  detail 
for  a  range  of  applications  by  various  authors.  Dede  (2009)  utilized  COM- 
SOL  Multiphysics  software  and  a  method  of  moving  asymptotes  optimizer 
to  Investigate  a  benchmark  heat  conduction  problem  Involving  internal  heat 
generation  and  a  heat  sink,  resulting  In  an  optimal  'branching'  structure. 
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The  author  then  extended  the  analysis  to  a  three-terminal  heat  transfer  and 
fluid  flow  device.  Chen  et  al.  (2010)  performed  multi-objective  topology  op¬ 
timization  for  finite  periodic  structures  and  considered  the  same  benchmark 
problem  In  combination  with  maximizing  the  stiffness  of  the  structure  under 
mechanical  loading.  Takezawa  et  al.  (2014)  considered  topology  optimization 
of  a  mechanical  structure  that  minimized  material  volume  under  both  strength 
and  thermal  conductivity  constraints.  Gersborg-Hansen  et  al.  (2006)  com¬ 
pared  heat  conduction  topology  optimization  results  modelled  by  the  finite 
element  method  and  the  finite  volume  method,  while  HasUnger  et  al.  (2002) 
used  the  homogenization  method  to  optimize  Isotropic  bi-materlal  conducting 
structures. 

One  particular  area  of  interest  has  been  that  of  multiple  heat  load  cases, 
where  a  number  of  heat  sources  act  on  the  structure  at  different  times,  lo¬ 
cations,  and  with  different  heat  boundary  conditions.  Zhuang  et  al.  (2007) 
used  a  level  set  method  to  explore  steady-state  heat  conduction  problems 
subject  to  multiple  heat  load  cases  with  the  design  objective  of  constructing 
an  effective  transport  path  for  heat  dissipation  under  a  given  volume  con¬ 
straint.  LI  et  al.  (1999)  Investigated  shape  and  topology  design  for  steady- 
state  heat  conduction  problems  using  the  ESO  method  for  both  single  and 
multiple  heat  load  cases.  Design-dependent  heat  loads,  or  loads  that  de¬ 
pend  on  the  material  distribution,  have  been  a  recent  topic  of  Investigation. 
Iga  et  al.  (2009)  considered  a  total  potential  energy  objective  to  determine 
optimally  conducting  structures  subject  to  design-dependent  boundary  condi¬ 
tions  of  heat  convection  and  Internal  heat  generation.  Gao  et  al.  (2008)  used 
a  modified  bidirectional  evolutionary  structural  optimization  (BESO)  method 
to  explore  steady-state  heat  conduction  under  both  design-independent  and 
design-dependent  heat  loads. 

Heat  conduction  also  plays  a  role  In  compliant  mechanism  topology  opti¬ 
mization,  where  the  mechanism  relies  on  the  devices  own  elastic  deformation 
to  transfer  a  motion  or  force.  Sigmund  (2001  a, b)  used  the  topology  optimiza¬ 
tion  method  to  design  thermally  and  electrothermally  driven  micro  actuators 
for  use  In  mlcroelectromechanlcal  systems  (MEMS).  In  these  systems  an  heat 
current  Is  converted  to  heat  which  causes  thermal  strain,  which  then  causes 
structural  deformation.  The  design  objective  Is  to  maximize  the  deforma¬ 
tion  of  a  workpiece  subject  to  a  constraint  on  material  volume,  heat  current, 
out  of  plane  displacement,  and  also  subject  to  the  heat,  thermal  and  elastic 
equilibrium  equations.  Yin  6  Ananthasuresh  (2002)  presented  a  design  pa¬ 
rameterization  scheme  for  topology  optimization  of  MEMS  made  of  multiple 
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Figure  1:  The  reference  domain 


materials  and  also  explored  design-dependent  boundary  conditions,  namely 
heat  convection  from  the  device  surfaces.  LI  et  al.  (2004)  designed  thermally 
actuated  compliant  mechanisms  that  consider  the  time-transient  effect  of  heat 
transfer  to  produce  the  localized  thermal  actuation. 


3  Optimization  of  thermal  conductivity 

We  consider  the  problem  of  steady-state  heat  conduction  with  no  convective 
heat  transfer.  Figure  1  depicts  a  body  defined  by  the  volume  Q  and  outer 
surface  F,  which  Is  subjected  to  a  prescribed  temperature  distribution  T*  on 
part  of  the  boundary  F t  and  the  normal  heat  flux  q*  Is  prescribed  on  the 
boundary  Vq. 


T=r  on  Vt 

q  •  n  =  q*  on  F^ 


(1) 

(2) 


The  governing  partial  differential  equation  for  steady-state  heat  conduction 
Is 

vMDV7)  +  S  =  0  (3) 

or  equivalently  In  Its  variational  form 


L 


VW  (DVr)  dO  + 


iySdO  =  0 


(4) 
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where  D  Is  the  conductivity  tensor,  S  Is  the  internal  heat  generation  and  W 
is  a  weight  function.  This  governing  equation,  also  known  as  Poisson  s  equa¬ 
tion,  may  be  used  to  describe  various  physics  problems  by  making  appropriate 
parameter  substitutions  (Donoso  6  Sigmund,  2004;  Huebner  et  al.,  2008). 

The  design  objective  of  maximizing  heat  conductance  is  achieved  by  min¬ 
imizing  the  thermal  compliance,  and  may  be  expressed  in  variational  form 
as 


min  : 

reT.p 

subject  to  : 


iiT] 


a[T,  W)  =  l(W). 

[  p(x)  c/0  <  V,, 
Jo 

p(x)e{o,i}. 


for  all  WgT. 


(5) 

(6) 


where 

a(T.W)^  [  VW  -lDVr)  dO  (7) 

Jo 

1{W)=  f  fVSdO-  f  Wq*dr  (8) 

Jo  Jr, 

and  T  denotes  the  space  of  kinematically  admissible  temperature  fields,  x 
Is  a  point  within  the  domain  0,  p(x)  Is  the  polntwlse  volume  fraction,  and 
Vf  Is  the  upper  bound  on  material  volume  fraction.  Using  a  finite  element 
discretization,  the  maximum  conductance  objective  may  be  expressed  as 


min  : 

X 

subject  to  : 


c(x)  =  T(x)TK(x)T(x) 
K(x)T(x)  =  F  =  F,  +  Fs 

0  <  X  <  1 


(9) 


where  x  is  the  vector  of  physical  element  densities,  c(x)  is  the  compliance, 
T(x)  is  the  global  temperature  vector,  K(x)  is  the  global  conductance  matrix, 
F  is  the  design  independent  global  thermal  load  vector  comprised  of  both  flux 
and  source  therms,  V'(x)  =  Yle=^  ^  fhe  material  volume,  Vq  is  the  design 
domain  volume,  Vf  is  the  prescribed  volume  fraction.  The  conductance  matrix 
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and  thermal  load  vector  are  assembled  from  their  element  contributions  as 
follows 


k,(x,)=  f  BjDe(x,)B,dO 

JCio 


(10) 


where  contains  the  shape  function  derivatives  and  De(xe)  Is  the  element 
conductivity  matrix  comprised  of  the  thermal  conductivity  values  as  a  function 
of  material  density.  For  an  Isotropic  material,  this  Is  given  by 

De(Xe)  =  de(Xe)\  (11) 

(12) 

where  de(Xe)  Is  the  element's  conductivity  as  a  function  of  element  density, 
and  I  Is  a  3x3  Identity  matrix.  Therefore,  the  element  conductance  matrix  may 
be  written  as 

kci'Xe)  =  de(Xe)  f  BjlB^dQ  (13) 

JOc 

=  de(Xe)ke  (14) 

where  Is  the  element  conductance  for  a  unit  thermal  conductivity.  The 

element  thermal  load  vector  Is  given  by 

r  =  f,^  +  f;  (15) 

=  N^’'q*dr+  f  N^’^SdQ  (16) 

-/r; 

Minimizing  the  compliance  minimizes  an  average  measure  of  temperature,  and 
therefore  maximizes  the  heat  heat  conduction. 

The  sensitivity  of  the  thermal  compliance  with  respect  to  physical  element 
densities  Xp  is 


ac(x)  ^ 


dXe 


(17) 


where  the  derivative  of  the  element  conductance  matrix  for  an  Isotropic  ma 
terlal  Is 

dk^(Xe)  dd{Xg)r 


dXe  dXe 

pXe  —  d^^^)  ke 

— _  d<°))  k,  for  RAMP 

L  (1  +  p(i  -  xe)y 


(18) 


for  SIMP, 
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where  p  Is  the  penalty  parameter  and  (1)  and  (0)  denote  the  two  material 
phases.  This  sensitivity  Implies  that  adding  conductive  material  to  the  design 
by  Increasing  the  element  density  will  decrease  the  compliance  and  therefore 
Improve  the  conductivity  of  the  structure.  This  Is  why  an  upper  limit  on 
volume  must  be  prescribed,  and  this  constraint  will  remain  active  throughout 
the  optimization  process.  As  with  the  previous  problem  formulations,  a  filtering 
step  converts  the  above  derivative  with  respect  to  physical  element  density 
to  a  derivative  with  respect  to  the  element  design  variable  Xp. 

We  have  tested  our  Implementation  with  the  aid  of  a  classical  heat  con¬ 
duction  problem  studied  by  numerous  authors  Including  Dede  (2009),  Chen 
et  al.  (2010)  and  Bendsoe  &  Sigmund  (2003).  The  problem  Involves  determin¬ 
ing  the  optimal  material  distribution  of  a  3D  structure  consisting  of  a  highly 
conductive  material  and  an  Insulator,  such  that  thermal  conductivity  Is  max¬ 
imized.  The  design  domain  undergoes  Internal  heat  generation,  with  a  heat 
sink  at  the  base  of  the  domain,  adiabatic  boundaries  and  a  constraint  on  the 
maximum  allowable  volume  of  the  conductive  material.  A  schematic  of  the 
problem  setup  Is  depleted  In  Figure  2. 


Adiabatic 
boundaries 
on  all  other 
faces 


Figure  2:  Boundary  conditions  for  the  heat  conduction  test  case 

Results  obtained  by  both  Dede  (2009)  and  Chen  et  al.  (2010)  are  shown  In 
Figure  3.  The  optimal  topology  exhibits  a  'branched'  structure  which  fans  out 
towards  the  boundaries  of  the  domain  so  that  heat  Is  drawn  from  the  domain 
down  to  the  heat  sink. 
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Figure  3:  Optimal  structure  obtained  by  (a)  Dede  (2009),  and  (b)  Chen  et  al. 

(2010) 


This  3D  heat  conduction  test  case  Is  Implemented  In  our  topology  op¬ 
timization  code  using  design  options  and  parameter  values  that  would  best 
replicate  the  problem  setup  used  by  previous  authors.  A  cubic  design  domain 
with  side  lengths  of  0.2  m  Is  discretized  Into  50x50x50  cubic  elements.  We 
only  analyze  one  quarter  of  the  domain  due  to  the  double  symmetry  of  the 
structure.  The  desired  volume  fraction  Is  set  to  =  0.30,  and  an  Initial 
homogeneous  distribution  of  material  Is  prescribed  where  the  element  den¬ 
sity  IS  set  to  the  given  volume  fraction.  Internal  heat  generation  Is  applied 
to  all  nodes  throughout  the  domain,  using  an  element  thermal  load  vector  of 
f/  =  0.01  *[1 ,1,1, 1,1, 1,1,1]^  W.  There  Is  s  heat  sink  In  a  rectangular  patch 
on  the  bottom  face  of  the  domain  where  the  temperature  is  set  to  T  =  0°C. 
For  this  problem  we  consider  two  materials,  both  of  which  exhibit  Isotropic 
thermal  conductivity.  The  highly  conductive  material  has  a  thermal  conduc¬ 
tivity  of  k  =  1  Wm“^°C“\  while  the  less  conductive  material  has  a  thermal 
conductivity  of  k  =  0.001  Wm“^°C“T  A  SIMP  Interpolation  scheme  is  used 
In  conjunction  with  a  continuation  scheme.  The  maximum  penalty  factor  Is 
set  to  pmax  =  3.  The  algorithm  uses  a  density  filter  with  a  filter  radius  of 
r  =  0.01  m  or  1.5  element  lengths.  The  method  of  moving  asymptotes  Is  the 
chosen  optimizer.  For  this  problem,  we  only  have  design  Independent  loading, 
therefore  the  sensitivity  Is  unconditionally  negative  and  the  MMA  Is  the  best 
optimizer  to  use.  The  GCMMA  method  will  yield  virtually  Identical  results 
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Figure  4:  Optimal  structure  obtained  using  our  topology  optimization  code 

for  this  particular  application,  albeit  with  a  greater  computational  time. 

Figure  4  shows  the  optimized  structure  obtained  using  our  topology  opti¬ 
mization  code.  It  appears  that  these  results  are  In  good  agreement  with  those 
presented  In  Figure  3.  Iteration  history  of  both  the  thermal  compliance  and 
volume  fraction  are  shown  In  Figure  5.  The  compliance  converges  by  approx¬ 
imately  50  Iterations  while  the  upper  volume  fraction  limit  of  vf  =  0.30  Is 
attained.  Possible  reasons  for  slight  discrepancies  between  the  three  opti¬ 
mized  structures  Include  use  of  different  material  properties  values,  thermal 
load,  filter,  filter  radius,  plotting  threshold,  continuation  scheme,  solver,  and 
optimization  algorithm. 

4  Multi-objective  optimization 

For  the  specific  problem  of  maximizing  heat  conduction  while  ensuring  struc¬ 
tural  Integrity  In  propellants,  promising  results  are  expected  by  combining 
the  maximum  structural  stiffness  and  the  maximum  conduction  objective  func¬ 
tions.  The  problem  formulation  Is  very  similar  to  that  presented  by  Chen  et 
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Figure  5:  Plot  showing  thermal  compliance  and  volume  fraction  vs  Iterations 


al.  (2010)  and  Is  shown  below: 


U^KsU 

mm  :  c(x)  = 


subject  to  : 


m 

Vo 


=  f 


Kc^  =  Fe 
KsU  =  Fs 


0  <  X  <  1 


(19) 


where  Wc  and  Ws  are  the  weighting  factors  for  the  maximum  conductivity 
and  maximum  stiffness  objectives  (Wc  -I-  =  1).  Kc  and  Kj  are  the  global 

conductivity  and  stiffness  matrices  respectively.  The  Individual  objectives  are 
normalized  by  their  maximum  values,  and 

Our  computational  methodology  and  Implementation  of  the  multi-objective 
case  was  verified  by  replicating  the  results  presented  by  Chen  et  al.  (2010). 
Chen  et  al.  developed  a  topology  optimization  algorithm  for  multifunctional 
3D  finite  periodic  structures,  simultaneously  addressing  the  maximum  stiff¬ 
ness  and  maximum  conductivity  criteria  using  a  weighted  average  method. 
The  system  setup  is  very  similar  to  our  model  where  the  design  objectives 
are  equivalent  while  the  loading  and  boundary  conditions  differ,  making  this 
example  an  excellent  test  case.  The  mathematical  formulation  Is  analogous 
to  Equation  19. 
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Applied 


Adiabatic 

boundaries 

Zero 

displacement 


Figure  6:  Loading  and  boundary  conditions 


The  following  parameters  are  used  In  the  optimization  algorithm 

•  Penalty  factor,  p  =  3 

•  Youngs  modulus,  E  =  1 

•  Poisson  s  ratio,  v  =  0.3 

•  Conductivity,  k  =  1 

•  Volume  fraction,  Vf  =  0.3 

•  Size  of  domain  =  80  x  80  x  80  elements 

Zero  displacement  Is  prescribed  at  the  bottom  four  corners  of  the  domain, 
while  a  unit  force  is  applied  In  an  upwards  direction  on  the  center  of  the  top 
surface.  The  design  domain  Is  heated  evenly  at  all  nodes  and  a  heat  sink  Is 
located  In  a  square  section  at  the  center  of  the  bottom  surface  (Figure  6). 

The  results  of  Chen  et  al.  are  provided  In  Figure  7.  The  Pareto  front 
shows  how  the  competing  design  objectives  influence  the  resulting  topologies 
as  the  weights  are  varied  from  =  0  (the  full  conductance  case)  to  W5  =  1 
(the  full  stiffness  case).  When  conduction  dominates,  the  optimal  design  Is  a 
doubly-symmetrlc  tree-llke  configuration  with  numerous  fine  twigs,  and  when 
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stiffness  dominates  we  obtain  a  structure  with  bars  connecting  the  fixed  points 
to  the  point  of  loading.  The  structural  topologies  obtained  using  our  code  are 
shown  In  Figure  8,  and  depict  a  similar  progression  of  structures  as  those 
obtained  by  Chen  et  al.  Slight  discrepancies  Include  the  structures  In  Figure 
8c  having  fewer  small  branches  and  differing  weight  values  for  visually  similar 
designs.  The  latter  can  be  attributed  to  the  objectives  being  normalized  by 
different  maximum  values,  while  other  causes  for  the  differences  Include  the 
use  of  a  different  filter,  filter  radius,  plotting  threshold,  continuation  scheme 
and  solver,  all  of  which  were  not  specified  In  the  paper.  Overall,  our  results 
appear  In  agreement  with  those  presented  In  Chen  et  al,  thereby  verifying 
our  computational  methodology  and  Implementation. 
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(a)  Ws  =  1  (b)  Ws  =  0 


Pareto  front 


Nonn  conduction  objective 


(c)  Pareto  front  and  corresponding  structural  topologies 


Figure  7:  Topology  optimization  results  produced  by  Chen  et  al.  (2010). 
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(a)  Ws  =  1 


(b)  ws  =  0 


Pareto  Front 
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Norm  conduction  objective 

(c)  Pareto  front  and  corresponding  structural  topologies 
Figure  8:  Topology  optimization  results  produced  using  our  algorithm 
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Archival  Publications  (published)  during  reporting  period: 

JOURNAL  PUBLICATIONS 

•  Theoretical  determination  of  anisotropic  thermal  conductivity  for  crystalline  1 ,3,5-triamino-2,4,6- 
trinitrobenzene  (TATB),  Matthew  P.  Kroonblawd  and  Thomas  D. 

Sewell,  J.  Chem.  Phys.  1 39,  074503  (2013). 

•  Theoretical  determination  of  anisotropic  thermal  conductivity  for  initially  defect-free  and  defective  TATB 
single  crystals,  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell, 

J.  Chem.  Phys.  141, 184501  (2014). 

•  Generalised  stacking  fault  energies  in  the  basal  plane  of  triclinic  molecular  crystal  1 ,3,5-triamino-2,4,6- 
rinitrobenzene  (TATB),  Nithin  Mathew  and  Thomas  D.  Sewell, 

Philosophical  Magazine  94,  424  (2015). 

•  Strategies  for  non-uniform  sampling  of  molecular  dynamics  phase  space  trajectories  of  relaxation 
phenomena,  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell,  Computer  Physics  Communications  196, 143 
(2015). 

•  Anisotropy  in  surface-initiated  melting  of  the  triclinic  molecular  crystal  1 ,3,5-triamino-2,4,6-trinitrobenzene: 

A  molecular  dynamics  study,  Nithin  Mathew,  Thomas  D.  Sewell,  and  Donald  L.  Thompson,  J.  Chem.  Phys. 
143,094706  (2015). 

•  Predicted  anisotropic  thermal  conductivity  for  crystalline  1 ,3,5-triamino-2,4,6-trinitobenzene  (TATB): 
Temperature  and  pressure  dependence  and  sensitivity  to 

intramolecular  force  field  terms,  Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell, Propellants,  Explosives, 
Pyrotechnics  41 , 502  (2016) 

•  Nanoindentation  of  the  triclinic  molecular  crystal  1 ,3,5-triamino-2,4,6-trinitrobenzene:  A  molecular 
dynamics  study,  Nithin  Mathew  and  Thomas  D.  Sewell,J.  Phys.  Chem.  C  1 20,  8266  (2016). 

•  Anisotropic  relaxation  of  idealized  hot  spots  in  crystalline  1 ,3,5-triamino-2,4,6-trinitrobenzene  (TATB), 
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Matthew  P.  Kroonblawd  and  Thomas  D.  Sewell,  J.  Phys. 

Chem.C  120, 17214  (2016). 

PAPERS  IN  PREPARATION  (meaning  a  complete  manuscript  is  being  finalized): 

•  A  density  functional  theory  study  of  Al  +  C02  reactions,  J.  D.  Veals,  R.  Chitsazi,  Y.  Shi,  T.  D.  Sewell,  and 
D.  L.  Thompson 

•  Predicted  melt  curve  and  liquid  state  transport  properties  of  TATB  from  molecular  dynamics  simulations, 
Nithin  Mathew,  Matthew  P.  Kroonblawd,  Thomas  D.  Sewell, 

and  Donald  L.  Thompson 

INVITED  PRESENTATIONS 

•  Matthew  P.  Kroonblawd,  Physical  and  Life  Sciences  Directorate,  Lawrence  Livermore  National 
Laboratory,  January  22,  2016 

•  Nithin  Mathew,  Theoretical  Division,  Los  Alamos  National  Laboratory,  June  21 , 2016 

•  Thomas  S.  Sewell,  Applied  Research  Institute,  University  of  Illinois  at  Urbana-Champaign,  July  14,  2014 

•  Thomas  D.  Sewell,  2016  Gordon  Conference  on  Energetic  Materials,  June  7,  2016 

•  Thomas  D.  Sewell,  Los  Alamos  National  Laboratory  MaRIE  ""  Workshop,  Probing  Dynamic  Processes  in 
Soft  Materials  Using  Advanced  Light  Sources,  July  27,  2016  (post  end-date,  but  covering  much  of  the  work 
accomplished  during  project  lifetime) 

•  Thomas  D.  Sewell,  2016  Lawrence  Livermore  National  Laboratory  "Mesoscale  Modeling  of  Explosives 
Initiation  Workshop,"  October  1 3,  201 6  (post  end-date,  but  covering  much  of  the  work  accomplished  during 
project  lifetime) 

•  Thomas  D.  Sewell,  To  be  presented  at  a  Los  Alamos  National  Laboratory  Workshop  "The  Behavior, 
Fabrication  and  Promise  of  Intentionally  Structured  Energetics",  to  be 

held  10-12  January  2017  (post  end-date,  but  covering  much  of  the  work  accomplished  during  project 
lifetime) 
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